Johannes Van Der Pol

Introduction to network modeling using Exponential Random Graph models (ERGM):

theory and an application using R-project

Introduction

Networks are representations of relational data. Nodes represent entities while the links connecting them represent any form of interaction or connection between the entities. A large diversity of networks exists ranging from networks of social contacts between individuals to inventor networks, collaboration networks, financial networks and so on. These networks attract the interest of researchers who want to explain the structure of these networks. In other words, the would like to know why certain agents are central in a network and why other are at the periphery, why some networks are densely interconnected and why others are sparse. In essence, we want to know why firms a is connected to firm b, or what it the probability that firm a will connect to firm c and how does this impact the overall structure of the network. This allows us to identify if technological proximity between firms has a significant impact on the choice of collaborator between firms and if this strategy has a significant impact on the structure of the network.

Financial support from IdEx Bordeaux and valuable comments and suggestions from Murat Yildizoglu are gratefully acknowledged.

The motivations for link creation cannot be observed directly from a visual representation of interactions, nor are they clear from glimpsing at a database containing relational data. In order to identify the motivations for entities to create links and identify the global network structure, a more in-depth analysis is required. Methods such as, Interpretive Structural Modeling, Total Interpretive Structural Modeling, and Graph Theory Matrix allow for an analysis of the interactions between variables but do not allow to for conclusions on the structure of the networks we aim at analysing. In other words, the methods allow for identifying a correlation between technological proximity and social proximity but do not explain the structure of the network that allows this conclusion. Econometric analysis could shed more light on the motivations behind an observed link through logistic regressions. The probability of a link could be explained by a number of variables. There is one important limitation to this method. Due to the hypothesis of independence of the observations the probability of a link between two nodes can never be explained by the presence of another link inside the network. It is feasible that a link between two nodes exists only because of the presence of other links in the network. Take for instance the idea that John and Mary are connected solely because they have a common contact: Paul. Methods such as Block-models cannot account for the impact of the structure of the network on the probability of a link. ERGM models are modified logistic regressions that allow for the probability of a link to depend upon the presence of other links inside the network (amongst other variables). ERGMs are able to take into account directed interactions as well as weighted interactions between nodes. The latter is an important point, since it allows an analysis beyond a simple binary relation between nodes approaching the idea of fuzzy logic. An ERGM identifies the probability distribution of a network so that it can generate large samples of networks. The samples are then used to infer on the odds of a particular link inside a network. Applications for this method are numerous in many fields of research as shown by the increasing trend in the number of publications using ERGM models (see figure 1). In economics the number of published papers appears to be relatively low when compared to the other social sciences. Only 19 published papers could be found in the Scopus database (and even less in the web of science database). The topics are however quite diverse: knowledge sharing in organizations [START_REF] Caimo | Knowledge Sharing in Organizations: A Bayesian Analysis of the Role of Reciprocity and Formal Structure[END_REF], GDP targeting [START_REF] Belongia | A working solution to the question of nominal GDP targeting // Macroeconomic Dynamics[END_REF], alliance networks [START_REF] Cranmer | Complex Dependencies in the Alliance Network // Conflict Management and Peace Science[END_REF][START_REF] Lomi | Relational collaboration among spatial multipoint competitors // Social Networks[END_REF][START_REF] Lomi | Networks in markets and the propensity of companies to collaborate: An empirical test of three mechanisms[END_REF][START_REF] Tom | Explaining the structure of inter-organizational networks using exponential random graph models // Industry and Innovation[END_REF] and geographic proximity [START_REF] Wal | The dynamics of the inventor network in German biotechnology: geographic proximity versus triadic closure[END_REF]. Since ERGMs allow for hypothesis testing, they can be put into use rather quickly within existing theoretical frameworks adding the possibility to analyse relational data in addition to normal data.

The growing interest, and development of a theory of economic networks, provides a fertile ground for the use of ERGM models from the geography of innovation to venture capital investments. The aim of this paper is to provide an overview of the basic statistical theory behind ERGM models, which will be dealt with in the first section. Section 2 discusses the concept of dependence and the explanatory variables that can be included in the models. Section 3 discusses estimation methods while section 4 provides the R scripts and the interpretation of an example using data for the French aerospace sector alliance network.

1 Theory

1.1 The canonical form of ERGM models The aim of an ERGM is to identify the processes that influence link creation. The researcher includes variables in the model that are hypothesised to explain the observed network, the ERGM will provide information relative to the statistical significance of the included variable much like a standard linear regression.

It is useful at this point to explain that sub-structures of a network can (and are predominantly) used as explanatory variables. Substructures are small graphs contained inside the network. Examples can be found in figure 2. The presence of some of these structures reflects certain link creation phenomena. A random network, i.e a network in which links are created at random, show a low number of triangles. A triangle is an interconnection of three nodes, the smallest possible complete subgraph. The presence of triangles in an empiric network bares witness that there is a process that generates triangles that is not the result of random link creation e.g, a tendency to create a link between common friends. A network with a small number of large stars and a large number of small stars can be the results of having a small number of very popular nodes. This is found in citation networks as well as lexicographical networks. Including sub-structures allows the modeling of certain processes as would any other variable.

In an ERGM we can find two types of explanatory variables: structural and node or edge-level variables. The latter come from other data sources and can be for example age, size of a firm, proximity, gender and so forth. The structural variables can contain indicators such as triadic closure, degree distribution and subgraphs.

The odds of a link

With this in mind the probability that one would observe a link between any two nodes i and j in a given network is proportional to a given set of explanatory variables:

p(G ij = 1) = θ 1 • X 1 + θ 2 • X 2 + . . . + θ n • X n (1)
We note G a graph and ij the focal nodes. G ij = 1 means that a link exists between nodes i and j in graph G, G ij = 0 implies the absence of a link between nodes i and j1 . θ is a vector of parameters and X a vector containing the variables. Equation 1 gives us the probability of a single link in graph G. Since nothing guarantees that the probability stays within [0, 1] the equation needs to be somewhat adjusted. We start by transforming the probability into an odds ratio:

odds(G ij = 1) = p(G ij = 1) 1 -p(G ij = 1) = p(G ij = 1) p(G ij = 0) (2)
Now, in equation 2 we notice that the odds of a link tend towards zero when the probability of a link tends towards one, while tending towards -∞ when the probability tends towards zero. A final modification is required to ensure that the probability remains between the bounds of 0 and 1, this is accomplished using the natural logarithm:

log(odds(G ij = 1)) = logit(p(G ij = 1)) = log p(G ij = 1) p(G ij = 0) (3)
The probability is now bounded by 0 and 1. If we suppose that the probability of a link is explained by a vector of n variables accompanied by their respective parameters (θ 1 . . . θ n ) then we can write:

logit(p(G ij = 1)) = θ • X (4)
So in our example we have:

logit(p(G ij = 1)) = θ 1 • X 1 + . . . θ n • X n (5) p(G ij = 1) = exp{θ 1 • X 1 + . . . θ n • X n } (6)
This gives us the logit which supplies the marginal probability of a tie between nodes i and j. Using these equations the probabilty of a tie between i and j would be independant from the probability between i and k. Indeed, a logistic regression works under the hypothesis of independence of observations. In the case of networks observations are not independent. For instance, common friends tend to connect more in social networks, common collaborators have a higher tendency towards collaboration. A model that aims at explaining a network structure should be able to include these tie formation processes. We hence modify the initial equations to include the network structure as observed before the link. This modification is introduced by Strauss, [START_REF] David | Pseudolikelihood estimation for social networks[END_REF]. We note G c ij the network without link ij:

odds(G ij = 1) = p(G ij = 1|G c ij ) 1 -p(G ij = 1|G c ij ) = p(G ij = 1|G c ij ) p(G ij = 0|G c ij ) (7)
In equation 7 the odds of a link between nodes i and j now depends on the structure of the network before a link between i and j is created (noted by |G c ij ). The probabilities are now conditional. We discussed previously that some of the variables in the model can be subgraphs. The manner in which these are included in the model is simply by the count of these sub-structures. In other words, the value of the variable triangles is the number of triangles in the network. The same is true for stars, circuits and shared partners. This has as a consequence that the counts of these variables are not the same when a link between two nodes is present or absent. For instance the number of edges changes by one. This means that we need to differentiate between the value of the variables when a link is present and when it is absent. We hence note the vector of variables v(G + ij ) when a link between i and j is added (hence the " + ") and v(G - ij ) when the link is absent. By including this differentiation we can rewrite equation 7 using the result in equation 6:

odds(G ij = 1) = p(G ij = 1|G c ij ) p(G ij = 0|G c ij ) = exp{θ • v(G + ij )} exp{θ • v(G - ij )} (8)
Where v(G + ij represents the vector of variables in the network with the link between i and j present and v(G - ij the vector of variables with no link between i and j. With some basic algebra we can develop the previous equation a bit further:

exp{θ • v(G + ij )} exp{θ • v(G - ij )} = exp{θ • v(G + ij )} • exp{-θ • v(G - ij )} (9) = exp{θ (v(G + ij -v(G - ij )} (10)
When developing the vector of variables we have:

= exp{θ 1 • (v 1 (G + ij ) -v 1 (G - ij )) + . . . + θ n • (v n (G + ij ) -v n (G - ij ))} (11)
Equation 11 shows that each parameter of the model is associated not with the counts of sub-structure but with the difference in counts. The difference from having an extra link, and the absence of said link. In essence

(v 1 (G + ij ) -v 1 (G - ij ))
represents the variation in the number of counts of network statistic 1 that result from the additional link. The variables are hence referred to as "change statistics". In order to remove the exponential from the right hand side of the equation we apply the logarithm:

log exp{θ • v(G + ij )} exp{θ • v(G - ij )} = θ 1 • (v 1 (G + ij ) -v 1 (G - ij )) + . . . + θ n • (v n (G + ij ) -v n (G - ij )) (12)
So we can rewrite equation 12 noting v 1 (∆ 1 G ij ) the change statistic for a link between i and j for variable 1 as follows:

log exp{θ • v(G + ij )} exp{θ • v(G - ij )} = θ 1 • v 1 (∆ 1 G ij ) + . . . + θ n • v n (∆ n G ij ) (13)
Each variable now accounts for the change in counts of network statistics. It is important to remind us that equation 13 accounts for the odds of one edge in the network while we are interested in the probability for the whole network. Following [START_REF] Julian | Nearest-neighbour systems and the auto-logistic model for binary data[END_REF] we can invoke here the Hammersley-Clifford theorem. Since this theorem is based on concepts out of the reach and the purpose of this document we will not detail the theorem. For a detailed explanation please refer to (Hammersley, Clifford, 1971). The theorem states that the probability of a network can be defined solely by the counts of subgraphs. This is important because it tells us that all we have to do is identify the correct subgraphs to ensure that a model of the network structure can be found. The more accurate the subgraphs to more reliable the inference of additional covariates.

The probability distribution of a network

The Hammersley-Clifford theorem states that the probability of a graph can be identified solely by counts of subgraphs. As such, we know that the probability is proportional to these variables. Since we have an observed network that we wish to replicate we look for the probability that the network generated by the model (X) is identical to the observed network (x). The logarithm is applied to bound the probability:

log(p(X = x)) ∝ θ • v(G) (14) p(X = x) ∝ exp{θ • v(G)} (15)
The right-hand side of the equation now needs to be normalized in order to obtain a proper probability distribution. The normalization is not straightforward, indeed, in order to normalize the probability of a network one needs to normalize by all possible networks with the same number of nodes:

p(X = x) = exp{θ • v(G)} y∈Y exp{θ • v(G)} (16)
With y a possible network structure in the set of all possible networks Y . The numerator is normalized by the sum of the parameters over all possible network structures. Note that this number is large. For a network with n nodes the number of possible graphs is 2 n(n-1) 2

. So even for a graph with 10 nodes there are 35184372088832 possible graphs. The major problem to overcome with ERGMs is exactly this normalizing constant.

With some simple algebra we find a general form for this model. Using θ as a vector of parameters and v(G) a vector of variables for network G:

p(X = x) = exp{θ • v(G)} exp{log( y∈Y exp{θ • v(G)}}) (17) p(X = x) = 1 ψ(θ) • exp{θ • v(G)} = exp{θ • v(G) -ψ(θ)} (18)
Equation 18 is the most general and commonly used form of the model [START_REF] Dean | Exponential random graph models for social networks: Theory, methods, and applications[END_REF]. Equation 18 also gives the canonical form of an ERGM model. Since the density of the random variable (the network structure) has the particular form in equation 18 it is referred to as an exponential family. In addition, since the structures are represented by a random variable, they are random graphs.

Putting both elements together and this results in an Exponential Family Random Graph. Since ERGM is easier to pronounce than EFRGM, the models are referred to as ERGM2 .

The canonical form gives the equation we wish to estimate. However, before we tackle the question of estimation we need to explore in more detail the variables that we would want to include in the model. We have stated previously that counts of subgraphs can be used as explanatory variables. The following section will explain which particular subgraphs are to be included in a model.

The dependence assumption

The previous section has shown that ERGMs are capable of providing conditional probabilities for links. This dependence assumption is important because it allows the researcher to study different phenomena that rule the formation of networks. This section will show how the hypothesis of dependence of links is connected to the choice of subgraphs that may be included in an analysis.

The Markovian hypothesis

Links between nodes rarely appear at random, agents have specific reasons to connect with one agent rather than another. The motivations behind interactions are numerous and complex and have been subject to scrutiny from researchers in different strands of the social sciences. This research shows that people or firms with common acquaintances or collaborators will have a higher tendency to cooperate for example. This makes sense for two reasons, first, having common partners means common working practices which have a positive impact on collaboration. Second, when searching for collaborations firms tend to rely on referrals. A collaboration already in place allows firms to observe in detail the efficiency of other firms, referrals that result from cooperation should hence be trustworthy. In addition cooperators of a firms have a higher probability to be in contact with each other since they are more likely to meet during social or professional events. The odds of a link depend on the neighborhood of the node and not on the entire rest of the graph. In more formal terms: two potential links X ij and X kl are in the same neighborhood if {i, j} ∩ {k, l} = 0. For instance, the dotted lines in Figure 3 represent potential links. The odds of this link will depend upon the nodes f and d have in common, in this case node a. This shows that not all subgraphs are compatible with analyzing this phenomenon. Any substructure that cannot account for a common node should not be included. Since it is hypothesized that only neighboring nodes impact the odds of a link, we seek node-level dependence. This level of dependence is also referred to as a nearestneighbors level of dependence3 or dyadic dependence [START_REF] Jenine | An introduction to exponential random graph modeling[END_REF]. For the purpose of the identification of relevant structures we need to find all structures that are subject to a Markov (or nearest-neighbor) level of dependence. When one identifies all the possible dependencies one can generate a dependence graph. The dependence graph for the complete Markov graph in our example can be found in figure 6. In red we find the dependence links for link a -b, it show the links on which a -b depends.

From this graph one can identify the substructures that comply with Markovian dependence. All subgraphs in the dependence graph can be included in a model to add Markovian dependence. For example, the dotted line between ab and ad represents a 2-star centered on agent a (figure 7).

Using the same method as for the 2-star one can also identify a triangle between the three agents on the links bd, bc and dc. The Markov model hence includes three configurations: edges, 2-stars and triangles. With the inclusion of these configurations the Markovian model takes the form:

In more complex graphs one could also identify 3-stars, 4-stars etc. It should be obvious here that the number of distinct 2-stars is large and it is near impossible to add a parameter for each distinct 2-star in the dependence graph. To reduce the number of variables a hypothesis is made that each type of configuration has the same probability of appearance, this allows for the inclusion of one parameter per substructure. In the case of Markov dependence the ERGM model would have the following form:

p(x = X|θ) = 1 ψ(θ) exp{θ E • v E (x) + θ S2 • v S2 + . . . + θ Sn-1 • v Sn-1 + θ ∆ • v ∆ } (19)
Where θ E is the parameter for the number of edges, θ S2 the parameter for the number of 2-stars and θ ∆ the parameter for the number of triangles. Note here that the model does not include simultaneously a 1-star and an edge parameter since they would be the same variable. With this model one is able to study if common nodes have a positive impact on the odds of link creation. In addition the combination of the 2-star parameter and the triangle parameter account for triadic closure effects. In other words, are triangles created because 3 nodes are connected at the same time, or are triangles formed by the closing of a 2-star.

Of course, Markovian dependence is only one of the possible levels of dependence. One can imagine higher levels of dependence, or even any level of dependence to be empirically relevant. One could suggest that firms evolving on the periphery of a network to have a higher probability to connect with firms in the center of the network than between them. The previous model was hence extended by [START_REF] Stanley | Logit models and logistic regressions for social networks: I. An introduction to Markov graphs andp[END_REF] to allow for a general level of conditional dependence giving the researcher a total liberty in the theories to test. Whatever the level of dependence chosen, the dependence graph gives the substructure that may be included [START_REF] Ove | Markov graphs[END_REF].

Higher levels of dependence

It is possible to assume that the Markovian level of dependence is not adequate or does not capture the full complexity of mechanisms of link creation. Links can be dependent without there being a common node involved. For instance, consider a case in which people work on the same floor in a company. The probability of a social link does not depend upon a common node but simply on the fact that they are geographically close, belong to the same community or have common cultural aspects [START_REF] Harrison | Identity and control: A structural theory of social action[END_REF]. In order to be able to model more complex aspects of social interactions and indeed even strategic interactions, one needs to be able to account for more structural aspects than stars and triangles (however potent in explanatory power these might be). The latter implies that the links are only dependent on each-other if nodes are part of a same neighborhood (neighborhood takes a broad definition here, it can be social, geographical or cultural). Due to the inclusion of general dependence the model is transformed to take the form:

p(x = X) = 1 ψ(θ) exp{ A∈M λ A • z A (x)} (20)
Where A indicates the neighborhood as part of the ensemble M of all possible neighborhoods. The parameter λ A will take a positive value when the probability of observing network x is increased. With the broad definition of "neighborhood" this model is able is almost limitless. The latter results in a problem, the model is too general.

We have seen in the previous subsection that the structures that can included in the ERGM model are defined by the dependency graph. In the case of a generalization of the dependency assumption, i.e all ties may be dependent upon all other ties, the dependency graph is a complete graph and all possible subgraphs can be considered as variables. This leaves a tremendous amount of parameters to be estimated. Of course, one can also include other types of variables to a model, such as age of the firm, geographic location, amount of public funds received etc. These variables are referred to as node variates or node attributes The addition of these attributes is introduced by Robins et al. (2001). The idea here is that links depend upon the attributes of the nodes they are connecting. In other words the probability that two nodes are connected depends upon the attribute of the node. These models are also called social selection models and take the following form:

p(Y = y | X = x) = 1 ψ(θ) • exp{ i θz(x) + θ a • z a (x, y)} (21)
Where the exponent a indicates parameters and structures relative to the inclusion of dyadic dependence for the purpose of social selection. In the same paper Robins et al. also described a social influence model in which the attributes of the nodes are a result of the structure of the network (nodes take a particular attribute according to the nodes in the neighborhood for example). In other words, the probability that a node variables takes a particular value depends upon the structure of the network and the values of this (or indeed any other) node-level variable.

We hence need to add a node variables to the model. Suppose we note Y i the value of a node variable for node i. This variable can be anything from a geographic region to the amount of public investment received by firms to the number of publications or patents. When this variable is included in the ERGM the model is written:

p(Y = y | X = x) = 1 ψ(θ i ) • exp{ i θz(y, x)} (22)
Just as it is possible to put values on nodes it is also possible to values on dyads. The question then is to know if the value of the dyad increases the probability of nodes being connected. Think of situations where we would like to know if the amount of investment between firms is related to cooperation or if technological proximity between firms induces collaboration. Note here that the difference between dyadic covariates and actor attributes resides in the value on the link between two nodes. In the case of proximity is refers to the proximity of both firms, it is hence not a firm-level variable. The value only makes sense when we consider firms two-by-two.

All the extensions made to the ERGM framework allow researchers to answer a large variety of questions about social and economic processes. Many other extensions which are beyond the scope of this document, but worth noting, are multivariate relations in networks [START_REF] Philippa | Logit models and logistic regressions for social networks: II[END_REF] and dynamic networks in which new ties depend upon the structure of the network at time t -1. It is also possible to model multivariate networks using ERGM. The idea is then that each link can exist in different matrixes, each corresponding to a different type of link (social, work, geography etc.). This extension allows researchers to study interplay between different networks and how each network is affected by the other networks.

Before looking at estimation methods for ERGM models one problem needs to be addressed: the degeneracy problem.

Solving degeneracy: Curved ERGMs

ERGM models are prone to degeneracy issues. When estimating the model the change statistics can behave in such a way that the large majority of the probability distribution is placed on either an empty or a full graph. As we will discuss in more detail a bit later, a simulation is performed to identify the probability distribution of a graph. This is done on a step by step basis, starting with an empty network and adding links one by one until a maximum likelihood is achieved. This probability distribution is a function of the change statistics and is thus impacted by the change in log-odds for an additional tie (for a given effect). In other words if an additional edge would create two new 2-stars, then the log-odds of observing that tie would increase by two multiplied by the parameter of the 2-star variable. A random graph model is considered stable if small changes in the parameter values result in small changes in the probabilistic structure of the model [START_REF] Mark | Assessing degeneracy in statistical models of social networks[END_REF]. When a new edge is added to the graph this not only increases the number of edges but might also increase the number of other sub-configurations that might be included in the model. The parameters of the model control the number of sub-graphs of a particular type that are expected in the graph. For instance a 2-star might be transformed into a triangle by the addition of an edge which also adds two 2-stars. A 2-star can become at 3-star and so on. This cascading effect will result in the simulation procedure jumping over the MLE and converge to a full graph. Lusher et al. (2012) [chapter 6] show that, in the case of the Markov (or triad model), the number of expected triangles increases as the parameter for triangles increases. They highlight a phase transition for certain values of the parameter where the number of expected triangles increases more than exponentially. This transition explain that the probability density distribution has high weights either on a (close to) empty graph or on a complete graph4 . This problem is increasingly present as the number of nodes increases. The larger the network the higher the number of substructures one can have. In order to avoid the model to put too much weight on the full graphs, [START_REF] Tom | New specifications for exponential random graph models // Sociological methodology[END_REF] propose to add several variables based on their concept of partial conditional dependence. The idea is to include a weighted degree distribution to the model, giving a high weight to low density while decreasing the weights as the degree increases. This reduces the impact of the high density variables responsible for the degeneracy of the initial model. Mathematically we can then write (using the notations of the initial paper):

u (d) α = n-1 k=0 e -αk d k (y) (23)
Where d k (y) is the number of nodes with degree k and α the parameter of the weights. This is referred to as the geometrically weighted degree distribution. The degree distribution can also be written as a function of the stars in the network. After all, a degree distribution is nothing more than a distribution of stars. Nodes with a degree of five are 5-stars, degree two are 2-stars and so forth. We can hence formulate the distribution as follows:

u s λ = S 2 - S 3 λ + S 4 λ 2 -. . . + (-1) n-2 • S n-1 λ n-3 = n-1 k=2 (-1) k • S k λ k-2 (24)
Where S k is a the number of stars of degree k and lambda the parameter. This method is referred to as Aternating k-stars. The difference between the geometrically weighted degree distribution and the K-stars is resides in the alternating signs. A large value of 3-stars is counterbalanced by a negative value for 4-stars due to the inverse sign of the parameter. The addition of the weights ensures that the change in change statistics stays small. Indeed [START_REF] Tom | New specifications for exponential random graph models // Sociological methodology[END_REF] show that the change statistics can be written:

z ij = -(1 -e -α )(e -α ỹ i+ + e -α ỹ j+ ) (25)
Where ỹ i+ is the density of firm i when the link between i and j is added. Equation 25shows that the value of the change statistic is reduced by the factor -(1 -e -α ). This factor hence ensures that the change statistics do not take too high values and result in a nested probability distribution. The inclusion of either the alternating k-stars or the geometrically weighted degrees transform the ERGM model into a Curved Exponential Random Graph Model [START_REF] Bradley | Defining the curvature of a statistical problem (with applications to second order efficiency)[END_REF].

All that is left to do now is estimate the model.

Estimation

Estimation allows for the identification of the parameters that maximize the likelihood of a graph. Since we only have one observation (the observed graph) a set of graphs from the same distribution is required. The set of graphs that may be generated by this procedure should have the observed graph as a central element to ensure a valid sample.

Markov Chain Monte Carlo

In the first section we identified the general form of an ERGM model (see equation 18). The odds of a graph were normalized by the sum of the parameters of all possible graphs. This leaves us with a constant to estimate which is near impossible. A workaround has to be found for ERGMs to be useful. A first development by [START_REF] Julian | Statistical analysis of non-lattice data // The statistician[END_REF]; Strauss, [START_REF] David | Pseudolikelihood estimation for social networks[END_REF] was to estimate the model using pseudo-likelihood estimation. The properties of this method are however not clear [START_REF] Tom | New specifications for exponential random graph models // Sociological methodology[END_REF][START_REF] Garry | An introduction to exponential random graph (p*) models for social networks // Social networks[END_REF] we shall hence focus here on more recent methods that are better understood.

A method for estimating the parameters of ERGMs using a sample is developed by [START_REF] Carolyn | A p* primer: Logit models for social networks // Social Networks[END_REF]; [START_REF] Tom | Markov chain Monte Carlo estimation of exponential random graph models[END_REF][START_REF] Charles | Constrained Monte Carlo maximum likelihood for dependent data[END_REF]. They estimate the model by Markov Chain Monte Carlo (MCMC) to find the maximum likelihood estimates (MLE). The idea is to extract a sample from a distribution that follows equation 18 asymptotically, not requiring the direct computation of the normalizing constant. Their paper points out that almost any maximum likelihood can be accomplished by a MCMC. A Markov chain is a sequence of random variables such that the value taken by the random variable only depends upon the value taken by the previous variable. We can hence consider a network in the form of an adjacency matrix in which each entry is a random variable. By switching the values of these variables to 0 or to 1 (adding or removing a link from the network) one can generate a sequence of graphs such that each graph only depends upon the previous graph. This would be a Markov chain. The hypothesis is then that if the value at step t is drawn from the correct distribution than so will the value at step t + 1. Unlike regular Monte-Carlo methods, the observations that are sampled are close to each-other since they vary by a single link. However, one would need a method for selecting which variable should change state in order to get closer to the MLE, this is done using the Metropolis-Hastings algorithm or the Gibbs sampler.

Metropolis-Hastings algorithm and the Gibbs sampler

The Metropolis-Hastings algorithm picks a variable at random and changes it's state. This results in either a new edge in the network or in the disappearance of an edge. The probability of the graph is then computed and only if the probability of the altered graph is higher than the previous one is the new graph retained for the next step. In other words the new graph is retained as long as the likelihood is increased:

min{1, p θ(x * ) p θ (x m-1 ) } (26)
This decision rule is called the Hastings ratio. The advantage of this ratio is that it does not include the normalizing constant ψ(θ).

Since the Markov chain starts at 0, a burn-in in needed to remove part of the chain to identify if the chain has converged or not (the burn-in can be parameterised in most software).

The steps taken by the Metropolis-Hastings algorithm are quite small. These small steps are implemented in order to avoid overstepping the global optimum which can easily happen in the case of larger parameter spaces. Other methods allow for bigger steps and as such converge faster and need a lower burn-in. The risk of larger steps is however overstepping the global optimum and convergence towards other local optima. The Metropolis-Hastings algorithm may be slower than others but is more precise in it's estimation.

Some programs use the Gibbs sampler, which is a special case of the Metropolishastings algorithm [START_REF] David | Inference in curved exponential family models for networks[END_REF]. The difference between Gibbs and Metropolis-Hastings resides in the chosen step. In the case of the Gibbs sampler, the state of each element in the vector of parameters is chosen and updated conditionally on the state of the other parameters. This means that if this decision rule was implemented in the Metropolis-Hastings algorithm the probability that the change is retained is equal to one. This makes the Gibbs sampler a relatively fast method. This sampling method is used by the different algorithms that are used to estimate ERGM models. Both methods allow for the generation of a sample of graphs that can be used for inference. The sample of graphs is obtained by varying not the parameters but the variables of the model until it is centered around the observed graph. Now that a sample of graphs has been found we need to estimate the parameters of the model. Two of the most widely used estimation algorithms, the "Stepping" and "Robbins-Monro" algorithm will now be reviewed.

The "Stepping algorithm"

This method introduced by Hummel et al. (2012). It has the advantage of approaching the MLE directly while "Robbins-Monro" does not. ERGM models are indeed estimated using the maximum likelihood method. Starting from the canonical ERGM form we define the log likelihood function as:

L(θ) = θ • v(G) -log(ψ(θ)) (27)
The problem here is the presence of the normalizing constant which cannot to computed. The improvements of this method over the previously one resides in the use of a log-normal approximation. The algorithm proposed here will converge towards the log-likelihood using a step-by-step method. The sampler used in with this estimation procedure is the metropolis-hastings sampler discussed previously. Once a sample of graphs has been identified the estimation algorithm is launched. Since the normalizing constant in equation 27 cannot be compute a workaround has to be found. The idea is to give starting parameters (θ 0 ).The log-likelihood ratio can then be written [START_REF] Hummel Ruth | Improving simulation-based algorithms for fitting ERGMs[END_REF]:

L(θ) -L(θ 0 ) = (θ -θ 0 ) T V (G) -logE θ0 [exp(θ -θ 0 ) T V (g)] (28)
Geyer, [START_REF] Charles | Constrained Monte Carlo maximum likelihood for dependent data[END_REF] point out that maximizing this ration by the means of a sample distribution of graphs generated with θ 0 only behaves well when θ is close to θ 0 . In other words one has to choose the correct starting point for the algorithm to find the MLE. The MLE solves the equation:

E θ v(G) = v(G obs ) (29)
The idea is to suppose that the MLE is not the observed value of the parameters but some point between the mean value parameterization and the observed value. A parameter γ defines the steps taken:

ωt = γ t • v(G) + (1 -γ t )ω (30)
Where ω t represents the estimate in the mean parameter space 5 .Ideally then, we would want γ = 1 so that the expected value of the parameters is the observed value. If this is the case the algorithm is considered to converge, this is shown in figure 8 which is the output of the R code. Once convergence is detected a large sample based on the parameters is computed and the MLE are estimated and gives as the final values.

Step 1 : Set the iteration number equal to 0 and choose initial parameters for vector η 0 .

Step 2 : Use MCMC to simulate graphs from the probability function for parameter vector η 0 .

Step 3 : Compute the mean of the sample.

Step 4 : Define a pseudo observation that is a convex combination of the mean of the sample and the observed value.

Step 5 : Replace the observed value by the pseudo observation.

1 I t e r a t i o n # 1 .

T r y i n g gamma= 0 . 1 7 2 I t e r a t i o n # 2 .

T r y i n g gamma= 0 . 1 4 3 I t e r a t i o n # 3 .

T r y i n g gamma= 0 . 1 6 4 I t e r a t i o n # 4 .

T r y i n g gamma= 0 . 2 2 5 I t e r a t i o n # 5 .

T r y i n g gamma= 0 . 2 5 6 I t e r a t i o n # 6 .

T r y i n g gamma= 0 . 3 4 7 I t e r a t i o n # 7 .

T r y i n g gamma= 0 . 3 1 8 I t e r a t i o n # 8 .

T r y i n g gamma= 0 . 4 6 9 I t e r a t i o n # 9 .

T r y i n g gamma= 0 . 7 4 10 I t e r a t i o n # 10 .

T r y i n g gamma= 0 . 9 7 11 I t e r a t i o n # 11 .

T r y i n g gamma= 1 12 I t e r a t i o n # 12 .

T r y i n g gamma= 1 13 Now e n d i n g w i t h one l a r g e s a m p l e f o r MLE. 14 E v a l u a t i n g l o g-l i k e l i h o o d a t t h e e s t i m a t e . U s i n g 20 b r i d g e s : 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 .

Stepping algorithm output

Robbins-Monro The Robbins-Monro algorithme is a stochastic approximation method introduced by Robbins, [START_REF] Herbert | A stochastic approximation method // The annals of mathematical statistics[END_REF] which is used by [START_REF] Tom | The statistical evaluation of social network dynamics // Sociological methodology[END_REF] and [START_REF] Tom | Markov chain Monte Carlo estimation of exponential random graph models[END_REF] to estimate ERGM models. Typically the method estimates:

E{Z θ } = 0 (31)
Where Z θ is a vector of parameters equal to u(Y ) -u 0 where u 0 is the observed value of the statistics. This allows us to rewrite the equation as a moment equation. The algorithm gives starting parameter values equal to the average of the parameter values. The initial parameters that will launch phase two are defined by [START_REF] Dean | Exponential random graph models for social networks: Theory, methods, and applications[END_REF]:

θ t+1 = θ t -a t • D -1 (z(x m ) -z(x obs )) ( 32
)
Where D is the co-variance matrix, the diagonal of this matrix will be used as the scaling matrix. a defines the convergence, it is set to a t = at-1 2 . The idea is that each step brings the values closer to the MLE. Hence large steps might result in the exceeding the MLE and divergence. The fact that the a r reduces in value with each step allows a smooth path to the MLE. As we move closer to the observed values of the statistics z(x m ) -z(x obs ) tends towards zero. The R output in figure 8 shows how the steps (a) start at a value of 0.1 and tend towards 0 with each iteration of the second phase of the algorithm. At the start of each step the starting parameters are considered to be the average values of the previous step. The number of iterations varies from model to model. The iterations stop once the trajectories of the generated statistics cross those of the observed ones [START_REF] Dean | Exponential random graph models for social networks: Theory, methods, and applications[END_REF]. The burn-in represents the number of simulations that are removed from the MCMC in order to make the chain "forget" the starting point. In other words it is to make sure the starting values do not impact the final result. Finally the algorithm checks for convergence using a convergence statistic. Just as in the case of the stepping algorithm one supposes that the MLE is reached when the distance between the observed values and the average of the simulated ones is close to 0. If there is no convergence than one can relaunch the estimation with as starting parameters the results of the previous simulation [START_REF] Dean | Exponential random graph models for social networks: Theory, methods, and applications[END_REF].

The largest difference between this method and the stepping method resides in two factors. First this method approaches an estimate of the MLE and does not evaluate the MLE function directly. Second, the steps are of a higher magnitude and can exceed the MLE if the starting values are close to the MLE. The use of either of the algorithms purely depends upon the model to be estimated. One algorithm might have better convergence in one case while the opposite can be true in another case. Note however that both use the Metropolis-Hastings method for the simulation of the MC. Now that we have discussed which variables can be included and how to estimate the parameters we will turn to an application using R.

1 Robbins-Monro a l g o r i t h m w i t h t h e t a _0 e q u a l t o : 2 e d g e s t r i a n g l e 3 -4.676219 1 . 4 5 6 3 8 0 4 P h a s e 1 : 13 i t e r a t i o n s ( i n t e r v a l =1 02 4) 5 P h a s e 1 c o m p l e t e ; e s t i m a t e d v a r i a n c e s a r e : 6 e d g e s t r i a n g l e 7 3 6 7 6 . 6 9 2 1 1 7 5 . 3 0 8 8 P h a s e 2 , s u b p h a s e 1 : a= 0 . 1 , 9 i t e r a t i o n s ( b u r n i n = 163 84) 9 t h e t a new : -4.66068075119643 t h e t a new : 1 . 3 9 9 3 2 8 1 3 3 9 0 9 5 4 20 P h a s e 3 : 20 i t e r a t i o n s ( i n t e r v a l =1 02 4) 21 E v a l u a t i n g l o g-l i k e l i h o o d a t t h e e s t i m a t e . We use here different R packages (Hunter et al., 2008a) , (Hunter et al., 2008b;[START_REF] Mark | statnet: Software Tools for Statistical Analysis of Network Data[END_REF]Butts, 2008).

The data (and hence the results of the estimations) are from the French aerospace collaboration network. Using patent data collaborations were identified which resulted in a network. The aim of the study is to understand if technological proximity played a significant role in the structuring of the collaboration network. We hence used a dyadic covariate called "proximity". The network contains 176 firms. The data used here was already in the form of an adjacency matrix and hence could be used directly. It is however also possible to use edgelists. Since the data needs to be transformed into a network object the network package will be needed. The latter is able to transform edgelists into adjacency matrices.

Since I'm using two measures of proximity I have two matrices and a matrix that includes the number of citations between firms. These need to be imported in the same manner as the network itself.

#We now t r a n s f o r m t h e i m p o r t e d d a t a i n t o n e t w o r k o b j e c t s w i t h t h e

p a c k a g e ' n e t w o r k ' 2 Network<-a s . m a t r i x ( Network ) 3 p r o x i m i t y<-a s . m a t r i x ( p r o x i m i t y ) 4 p r o x i m i t y . e<-a s . m a t r i x ( p r o x i m i t y . e ) We now have different objects to work with: the network and dyadic covariates in the form of networks. Note that the model constructed here is for un undirected network as shown by the option directed = F ALSE in the as.network function. We now have to decide which variables to include in the model. Let's start with a simple model containing only edges. We invoke here the ergm() function from the ergm package:

1 model<-ergm ( Network ~e d g e s )
This gives us the most basic form of an ERGM model, the estimation method defaults to Monte Carlo MLE [START_REF] Charles | Constrained Monte Carlo maximum likelihood for dependent data[END_REF]. The parameter for the variable edges has an estimated value of -3.8528. This means that the probability of two ties connecting is:

p(i → j) = exp(-3.85) 1 -exp(-3.85) = 0.02174241 (33)
Recall equation 11, this equation stated that the variables were change statistics. The parameter should hence be multiplied by the change in the number of subgraphs. In other words, if an additional edge creates three triangles the parameter should be multiplied by three. Since an additional edge only creates one new edge we do not multiply. Let's try the same but with only triangles as explanatory variable. The parameter for the variable triangles has an estimated value of -1.9146. This means that the log-odds of two nodes connecting is:

-1.9146 * δ triangles (34)
Where δ triangles gives the change in the number of triangles. Hence the log-odds depend upon the number of triangles that will be created by an additional tie: > If the link creates 1 triangle, the log-odds are 1 * -1.9146. The probability is then 0.1284 > If the link creates 2 triangles, the log-odds are 2 * -1.9146. The probability is then 0.0213

Let's see how we interpret estimates when we have more than one variable: the Markov model. In this case an additional edge, if it creates x 2-stars and y triangles, has log-odds:

-28.30 + x • 7.68 + y • 57.099

Since the model includes lower and higher order subgraphs (2-stars are a substructure of triangles) we can conclude here that triadic closure is significant in the network. In other words, 2-stars are closed to form triangles. Note however the values of the information criteria AIC and BIC, stating these values are high is an understatement. The model must suffer from some kind of problem. Subsection 4.3 will show how to deal with degeneracy and other problems. Now we want to see how proximity (technological in the case of this study) influences the structuring of the network. Let's start with a simple model using the edge covariate "proximity". The comment for adding an edge covariate is simple edgecov. Similarly node covariates (age, gender, type of firm, country...) are added using the nodecov() or nodef actor command. The node factor command is particularly useful since it allows to compare log-odds to a reference point. For example, one could categorise the R&D expenses of firms into low, average and high. The ERGM would then give the log-odd of 2 categories as compared to a third. In other words is a link more likely for firm with average expenses than for low ? This feature seems not to be available for edge covariates however.

The following model only contains edges and the edge covariate. A first point we can notice is that the AIC and BIC criteria are lower with the addition of the edgecovariate. The model is hence improved with the addition of the proximity parameter. Firms with proximity are 2.13 times more likely to connect in this network (citeris paribus)6 .

The probability of an additional edge is then positively impacted by the technological proximity. More specifically the average degree of the network is positively impacted by technological proximity. A network analysis performed on this network showed that the network has a scalefree structure. This information can be helpful in the modeling of the ERGM as we have information on the distribution of the degrees. The same is valid for any other information about the network, small-world properties, level of clustering or centrality distribution. The information provided by SNA allows a first understanding of the structural properties of the network that will allow for a more robust model once edge and nodal covariates are added. So if the network ha a scale-free structure the structure should be explained by the degree distribution. To check this we can add different degrees to the model. This can be done by using the command degree(). On can simple add one statistic for one degree, i.e degree(3), for the impact of nodes with degree 3, or add multiple degrees as was done in the following model. Note that the addition of multiple degrees is achieved by writing degree(a : b) to add all the degrees between a and b. The addition of these variables to the model once again decreases the AIC and the BIC, the model is hence enhanced. The structure of the network can be explained by a degree distribution. These models seem to work qui nicely. We discussed in the previous sections that models were prone to degeneracy and the solutions to this problem. Let's have a look at how we can model ERGMs with alternating k-stars and a geometrically weighted degree distribution.

Curved Exponential Random Graph Models

Studies show that the addition of weights on the degree distribution helps to avoid bi-modal distributions in the parameter space, i.e avoids the generated networks from being either full or close to empty. Different forms can be added to the R code. Since we have here an undirected network we can use either the alternating k-stars altkstar() or the geometrically weighted degree distribution gwdegree. For directed graph there are additional commands which work in similar manner as what we show here. We include here a statistic for the gwdegree with the option f ixed = T RU E. The latter means that we do not make an estimation of the scaling parameter, we want it to be equal to 1. The resulting model is hence not a curved ERGM. In order to have a curved exponential random graph model, the parameter that defines that we fixes in the previous code has to be estimates as well. In the following code we estimated the model with an edgewise shared partners variable. This variable is used to check for transitivity. By switching the option f ixed = T RU E to f ixed = F ALSE the model becomes curved. The results now include an estimate for the parameter alpha of the model. Note here that the parameter can only be interpreted if the gwesp statistic is significant. The interpretation of these estimates are much more complex than previously. The parameters need to be exponentiated to find λ that we saw in the equations [START_REF] David | Curved exponential family models for social networks // Social networks[END_REF]. We hence find e 0.88784 = 2.42. Since the parameter is positive we can conclude that transitivity is present.

Goodness of fit diagnostics

In order to check if a model is a good fit we use the mcmc.diagnostics command. This gives us a number of outputs, notable the matrix of correlations and p-values for both the individual parameters and the model as a whole. The p-values are high for the parameters and the model, we can hence conclude that the model is globally significant. In order to go into a bit more detail when it comes to the estimates, the commands also provides us with plots, see figures 10 and 11. We stated that an ERGM should fit the observed network perfectly, on average. This means that from the simulated networks we expect the average values to be those of the observed network. If this is not the case then the sample the model is based on does not come from the same distribution as the observed network. Figure 10 shows an example of what we want to observe. in the first graphs the values oscillate around the mean which is what we want. The graphs on the right hand show a centered distribution of the values, we hence conclude that the model is a good fit. A bad example can be found in figure 11. The graphs show that the distribution is not at all centered, and there is no oscillation around the mean. This model is hence a bad fit.

One can also study the goodness of fit using the gof () command. Using a plot command this provides a box plot (see figure12). The gof command can receive different parameters, one can chose to plot any number of variables and decide to increase or decrease the number of simulations to refine the results. The following code provides the GOF of the whole model (GOF = M odel) using 20 simulations (nsim = 10).

1 g o f _ model<-g o f ( model , GOF=~Model , nsim =2 0 )
The results give us a box plot per variable and a black line representing the observations on the empirical network. Since we want the mean of the simulations to be equal to the observed network, the dark line should coincide with the center of the boxplots (the vertical line in the boxplot representing the median of the distribution). This is the case here, the model is hence a good fit.

Improve bad models

The fitting of an ERGM is a trial and error procedure. If a model behaves badly there are a couple of parameters to change in order to improve results. Of course one should Once this is done and degeneracy is still observed one can start by switching estimation methods. One method might work better in one case than the other.

The estimation algorithm can be chosen in the control arguments of the ergm() command.

1 model<-ergm ( Network ~e d g e s , c o n t r o l = c o n t r o l . ergm ( main . method = " S t e p p i n g " ) ) 2 model<-ergm ( Network ~e d g e s , c o n t r o l = c o n t r o l . ergm ( main . method = " S t e p p i n g " ,MCMC. s a m p l e s i z e =70000 , MCMC. i n t e r v a l =5 00 0) )

The burn-in can also solve problems, the burin represents the number of iterations that are removed from the simulation. It other words, the higher the burnin the more the procedure forgets about it's initial parameters. Increasing this value hence allows for keeping only the latest values which might represent the real values better. This can be achieved by adding an option to the ergm. A second method of improving estimation would be to increase the sample size by changing the parameter M CM C.samplesize. This increase will result in having more precise estimates by an increase in the number of statistics drawn from the sample. This, of course, increases computation time. The ergm package includes multicore features that can help reduce computation time drastically. All this requires is the addition of some parameters to the control.ergm argument. Adding parallel = 4 notifies the package that the computer has 4 cores, parallel.type sets the type of multicore. For a regular computer this should be fixed to "P SOCK". This will distribute the computations over the 4 cores of the computer and hence increase speed.

1 a<-ergm ( Network ~e d g e c o v ( p r o x i m i t y 2 ) + t r i a n g l e s + a l t k s t a r ( 1 . 8 1 2 , f i x e d = FALSE ) , c o n t r o l = c o n t r o l . ergm ( main . method =c ( " S t e p p i n g " ) , p a r a l l e l =4 , p a r a l l e l . t y p e = "PSOCK" , MCMC. s a m p l e s i z e =20 00 0) )

Conclusion

With the interest for network analysis growing in economics, ERGM are an efficient tool for an in-depth analysis while checking for statistical significance. The complexity of social and economic phenomena are difficult to assess and even though existing methods such as block-models and classic logistic regressions do not allow an assessment of all the intricate rules that are at play. By taking into account the impact of the current structure of the network ERGM distinguish themselves from other methods of network analysis and allow for more accurate hypothesis testing. Being able to add relational data to an analyses allows for increasing the precision of existing models as well as the testing of new hypotheses.

Even though the models are powerful, the tools used for their analysis still need to be improved. The fast expansion of research on the subject of ERGMs is tackling issues improving the quality of ERGMs overall. One important issue is worth mentioning here, which is the size of the network under analysis. In the present paper the network had a size of 176 nodes, large networks can cause significant issues with computation time or the program simply refusing to run. Recent research is aiming at tackling this issue (Schmid, Desmarais, 2017;Bouranis et al., 2017). Many questions in economics are related to the strategies of link creation which are difficult to answer with other methods of analysis. Whether they are used to analyse social networks, collaboration networks, trade networks or financial networks, ERGMs can provide vital insights into the understanding of network dynamics. I should however point out at this point that the hypothesis of modeling networks with an exponential law should not be forgotten. Certain interactions might not be well suited for exponential laws, a notable exemple is citations networks. Research has shown that citation networks follow a power-law rather than an exponential law. For networks with an underlying power-law one should look at Dynamic Process Modeling which is based on a power-law.

Figure 1 :

 1 Figure 1: Evolution of the number of publications involving ERG models for all disciplines (statistics included) (source: Scopus)

  Figure 2: Network G

Figure 3 :

 3 Figure 3: Node level dependence illustration: the Markovian neighborhood

  Figure 4: Markov graph

Figure 8 :

 8 Figure 8: R output for the Stepping algorithm

  a s e 2 , s u b p h a s e 2 : a= 0 . 0 5 , 23 i t e r a t i o n s ( b u r n i n = 163 84) 12 t h e t a new : -4a s e 2 , s u b p h a s e 3 : a= 0 . 0 2 5 , 58 i t e r a t i o n s ( b u r n i n = 163 84) 15 t h e t a new : -4.62881856406474 16 t h e t a new : 1 . 4 1 4 0 5 5 9 3 9 6 6 0 4 2 17 P h a s e 2 , s u b p h a s e 4 : a= 0 . 0 1 2 5 , 146 i t e r a t i o n s ( b u r n i n = 163 84) 18 t h e t a new : -4.60985096388914 19

Figure 9 :

 9 Figure 9: R output for the Robbins-Monro algorithm

  1 # I m p o r t d a t a . 2 Network<-r e a d . t a b l e ( "ADJ_MATRIX . c s v " , s e p = " ; " )

  r a n s f o r m t o n e t w o r k f o r m a t 7 Network<-a s . n e t w o r k ( Network , d i r e c t e d = FALSE ) 8 p r o x i m i t y<-a s . n e t w o r k ( p r o x i m i t y , d i r e c t e d =FALSE ) 9 exp _ p r o x i m i t y<-a s . n e t w o r k ( exp _ p r o x i m i t y , d i r e c t e d =FALSE ) 10 c i t a t i o n <-a s . n e t w o r k ( c i t a t i o n , d i r e c t e d =TRUE)

  i m a t e S t d . E r r o r MCMC \%p-v a l u g n i f . c o d e s : 0 ' * * * ' 0 . 0 0 1 ' * * ' 0 . 0 1 ' * ' 0 . 0 5 ' . ' 0 . 1 l D e v i a n c e : 21349 on 15400 d e g r e e s o f f r e e d o m 16 R e s i d u a l D e v i a n c e : 3113 on 15399 d e g r e e s o f f r e e d a l l e r i s b e t t e r . )

  g n i f . c o d e s : 0 ' * * * ' 0 . 0 0 1 ' * * ' 0 . 0 1 ' * ' 0 . 0 5 ' . ' 0 . 1 ' ' 1 14 15 N u l l D e v i a n c e : 21349 on 15400 d e g r e e s o f f r e e d o m 16 R e s i d u a l D e v i a n c e : 8935 on 15399 d e g r e e s o f f r e e d o m 17 18 AIC : 8937 BIC : 8945 ( S m a l l e r i s b e t t e r . )

  04 * * * 14 ---15 S i g n i f . c o d e s : 0 ' * * ' 0 . 0 0 1 ' * * ' 0 . 0 1 ' * ' 0 . 0 5 ' . ' 0 . 1 ' ' 1 16 17 N u l l D e v i a n c e : 21349 on 15400 d e g r e e s o f f r e e d o m 18 R e s i d u a l D e v i a n c e : 35811780 on 15397 d e g r e e s o f f r e e d o a l l e r i s b e t t e r . )

  14 S i g n i f . c o d e s : 0 ' * * * ' 0 . 0 0 1 ' * * ' 0 . 0 1 ' * ' 0 . 0 5 ' . ' 0 . 1 ' ' 1 15 16 N u l l D e v i a n c e : 21349 on 15400 d e g r e e s o f f r e e d o m 17 R e s i d u a l D e v i a n c e : 3096 on 15398 d e g r e e s o f f r e e d o a l l e r i s b e t t e r . )

  17 e d g e c o v . p r o x i m i t y 0 l D e v i a n c e : 21349 on 15400 d e g r e e s o f f r e e d o m 22 R e s i d u a l D e v i a n c e : 2989 on 15393 d e g r e e s o f f r e e d a l l e r i s b e t t e r . )

  r m u l a : Network ~e d g e s + d e g r e e ( 2 : 6 ) + e d g e c o v ( p r o x i m i t y 2 ) 6 7 I t e r a t i o n s : 7 o u t o f 20 8 9 Monte C a r l o MLE R e s u l t s : 10 E s t i m a t e S t d . E r r o r MCMC \% p-v a l u g n i f . c o d e s : 0 ' * * * ' 0 . 0 0 1 ' * * ' 0 . 0 1 ' * ' 0 . 0 5 ' . ' 0 . 1 l D e v i a n c e : 21349 on 15400 d e g r e e s o f f r e e d o m 22 R e s i d u a l D e v i a n c e : 2989 on 15393 d e g r e e s o f f r e e d a l l e r i s b e t t e r . )

  r m u l a : Network ~e d g e s + t r i a n g l e s + e d g e c o v ( p r o x i m i t y 2 ) + g w d e g r * * * 16 ---17 S i g n i f . c o d e s : 0 ' * * * ' 0 . 0 0 1 ' * * ' 0 . 0 1 ' * ' 0 . 0 5 ' . ' 0 . 1 ' ' 1 18 19 N u l l D e v i a n c e : 21349 on 15400 d e g r e e s o f f r e e d o m 20 R e s i d u a l D e v i a n c e : 4170 on 15396 d e g r e e s o f f r e e d o a l l e r i s b e t t e r . )

  r m u l a : Network ~e d g e s + e d g e c o v ( p r o x i m i t y ) + gwesp ( a l p h a = 1 , f i x e d = FALSE ) 6 7 I t e r a t i o n s : NA 8 9 S t e p p i n g MLE R e s u l t s : 10 E s t i m a t e S t d . E r r o r MCMC % p-v a l l D e v i a n c e : 21349 on 15400 d e g r e e s o f f r e e d o m 19 R e s i d u a l D e v i a n c e : 2847 on 15396 d e g r e e s o f f r e e d

1

  Sample s t a t i s t i c s c r o s s -c o r r e l a t i o n s d i v i d u a l P-v a l u e s ( l o w e r = w o r s e ) : 7 k s t a r 2 e d g e c o v . p r o x i m i t i n t P-v a l u e ( l o w e r = w o r s e ) : 0 . 6 9 9 4 2 3 3 .

  Figure 10: Goodness of Fit diagnostics

Figure 11 :

 11 Figure 11: Goodness of Fit diagnostics, bad example

  Figure 12: GOF: boxplot analysis

Figure 13

 13 Figure 13: Triads

Figure 16 :

 16 Figure 14: K-stars

the values 0 and 1 refer to values found in an adjacency matrix, 1 indicating the presence of a link, 0 the absence

This type of model is also referred to as the P * family of models[START_REF] Carolyn | A p* primer: Logit models for social networks // Social Networks[END_REF][START_REF] Dean | Exponential random graph models for social networks: Theory, methods, and applications[END_REF].

Nearest-neighbor systems have been studied by[START_REF] Julian | Nearest-neighbour systems and the auto-logistic model for binary data[END_REF].

In addition, the degeneracy of the model can result in problems with the estimation procedures.

In the algorithm the initial values are chosen to be the MPLE

# I m p o r t t h e d y a d i c c o v a r i a t e s

p r o x i m i t y<-r e a d . t a b l e ( " P r o x i m i t y _ m a t r i x . c s v " , s e p = " ; " )

p r o x i m i t y . e<-r e a d . t a b l e ( " P r o x i m i t y _ m a t r i x _ exp . c s v " , s e p = " ; " , d e c = " , " )

c i t a t i o n <-r e a d . t a b l e ( " C i t a t i o n _ m a t r i x . c s v " , s e p = " ; " )

e 0.7575 since this is the odds and not the log-odds


Henri Buscail
email: henri.buscail@udamail.fr

Raphaël Rolland

Frédéric Riffard
email: frederic.riffard@udamail.fr

Christophe Issartel
email: christophe.issartel@udamail.fr

Cécile Combe
email: cecile.combe@cetim.fr

Pierre-François Cardey
email: pierre-francois.cardey@cetim.fr

Publication in the journal 'Oxidation of Metals' (Peer-reviewed) Influence of Sodium Hydroxide on the Austenitic 330Cb Alloy Oxidation at High Temperature

Keywords: NaOH, 330CB Alloy, SiO 2, Cristobalite

Sodium hydroxide is found on metallic conveyor belts as residue of the degreasing process of metallic screws put inside heat treatment furnaces at 900 °C. The 330Cb (Fe-34Ni-23Cr-1Nb-1.55Si) alloy is used to build these conveyor belts. With time, the amount of sodium hydroxide increases on the 330Cb steel surface. The aim of this work is to show the influence of the NaOH amount on the 330Cb oxidation during 48 h, at 900 °C. Results show that, on this austenitic steel, a protective silica scale formation is promoted by low oxygen containing gaseous environments and the high alloy silicon content. With low NaOH deposits the oxide layer is adherent and a SiO 2 cristobalite subscale is formed. It is also demonstrated that sodium combines with niobium to form NaNbO 3. When the NaOH amount exceeds 0.265 mg.cm -2 on the alloy surface, severe oxidation is observed leading to a fast growing FeCr 2 O 4 oxide scale.

INTRODUCTION

The 330Cb alloy (Fe-34Ni-23Cr-1Nb-1.55Si) is generally used in high temperature carburising and nitriding conditions. These environments induce severe degradation of the conveyor unit due to a rapid weakening of the alloy. To improve the alloy carburization resistance it is assumed that the presence of an adherent oxide scale acts as a carbon and nitrogen diffusion barrier. A previous study has shown that in the 800 -1000 °C temperature range, the 330Cb oxidation leads to a chromia scale acting as a good diffusion barrier under isothermal conditions [1]. The external manganese chromite subscale can limit the chromia scale evaporation at temperatures higher than 1000 °C [2]. Nevertheless, after cooling to room temperature, important oxide scale spallation occurred [1]. Gleeson demonstrated the significant effects of alloy composition on long-term, cyclicoxidation resistance. Each of the alloys exhibited scale spallation; however, the manner by which spallation occurred varied between the alloys [3]. In the literature, a Fe-35Ni-18Cr-2Si alloy has been studied at high temperature [4]. It has been shown that a Cr 2 O 3 scale is protective at high temperature. Nevertheless, Stevens has indicated that at temperatures higher than 1100 °C the chromium oxide CrO 3 evaporation occurs. Then, chromium depletion is observed in the oxide scale and in the alloy. On the basis of changes in the chemical composition of the oxide film, together with spalling, the failure of heating-element wire by the development of hot spots can be explained. At failure both long-term (low temperature < 1000 °C) and shorter-term tests (temperature > 1000 °C) were found to have oxidised to the same depth. This can be explained on the basis of a depletion of chromium, which occurred at the oxide-alloy interface, and becomes insufficient for a healing process. Gleeson's work stated that the HR-160 alloy exhibited complete spallation owing largely to its relatively high silicon content (2.75 wt.%) [3]. However, the silicon was also beneficial in promoting protective scale formation when the exposed alloy was subsequently oxidised. The HR-120 alloy showed the poorest cyclic-oxidation resistance, due to poor scale adhesion and the tendency of iron (33 wt.%) to oxidise. All of the alloys underwent internal oxidation and void formation. In most cases, the extent of internal attack was significantly greater than that of metal loss. Former studies have shown that the oxidation of the 330Cb alloy at 900 °C in air lead to a poorly adherent chromia scale and did not permit a continuous silica scale formation even though the alloy contains about 2 wt.% silicon [5][6][7]. Then, the protection against corrosive environments is not ensured. One way to act on the oxide scale formation consists of modifying the gaseous environment, using inert gases [8]. Previous works have shown the effect of sodium salts coating on the high-temperature oxidation of alloys. Sodium chloride and sodium sulphate induce accelerated oxidation of chromium and hot corrosion processes [9][10]. Sodium carbonate coated Nickel-base alloys oxidised in air at 900 °C also show fluxing reactions but the morphology of the scales show less perturbation compared to sodium sulphates, sodium chlorides or sodium nitrates [11].

In this work, chlorine is not present and we consider the presence of sodium hydroxide on metallic conveyor belts as residue of the degreasing process of metallic pieces put inside heat treatment furnaces at 900 °C. The 330Cb alloy is often used to build these conveyor belts. With time, the amount of sodium hydroxide can increase on the 330Cb alloy surface due to residue accumulation. Based on what is known about the oxidation behaviour of this alloy, this work will focus on the influence of sodium hydroxide coatings on the 330Cb oxidation during 48 h, at 900 °C. The low oxygen potential gaseous environment used, N 2 -5vol.% H 2 , simulates industrial conditions in the heat treatment furnaces. In order to explain the influence of NaOH on the 330Cb oxidation, we will examine the effect of various NaOH coatings on the 330Cb alloy oxidation at 900 °C.

EXPERIMENTAL PROCEDURES

The 330Cb alloy is an austenitic stainless steel. Its composition, in weight %, is given in Table 1. 1.5 mm thick cylindrical specimens of 12 mm diameter were abraded up to the 320-grit SiC paper, then degreased with ethanol and finally dried. All oxidation tests were performed at 900 °C. High temperature oxidation was performed during 48 h in flowing (8 l/h) nitrogen containing 5 vol.% hydrogen (N 2 -5vol.% H 2 ). The residual pO 2 = 15 ppmv contained in the gas was measured by an oxygen analyser (Elcowa GPR 1200MS). The oxide characterisation was realised by X-ray diffraction (XRD). XRD patterns were obtained by use of a Philips X'pert MPD diffractometer (copper radiation, λ α k = 0.15406 nm).

The XRD conditions were 2Θ scan, step 0.05° ranging from 10 to 80°, 8 s counting time. The oxide scale surface and cross-section morphologies have been observed in a JEOL 7600 scanning electron microscope (SEM) coupled with a LINK energy dispersive X-ray spectroscopy (EDXS). The EDXS point analyses were performed with an electron probe focused to a 1 µm spot. The sodium-hydroxide coatings were realised by using: 0.01; 0.1 and 0.3 mol.l -1 solutions. Immediately after dipping in the solution, the specimens were dried in warm flowing air at 40°C during 5 min. After coating, the specimens were weighted in order to determine the deposit weight per unit area. It corresponds respectively to: 7.10 -3 ; 6.8 10 -2 ; 0.262 mg.cm -2 . After drying, the coating was characterized by XRD and it was found that sodium hydroxide reacts with ambient carbon dioxide present in the air to form sodium carbonate Na 2 CO 3 (ICDD 18-1208). During oxidation, Na 2 CO 3 decomposes at 856 °C to form Na 2 O and this compound is reacting with the alloy at 900 °C. XRD patterns obtained on blank specimens before oxidation only show the presence of the autenitic structure (ICDD 03-1209). XRD results also indicate that no niobium containing intermetallics such as Nb 3 Ni 2 Si (ICDD 17-0097) could be detected on the alloy surface before oxidation.

RESULTS

On Figure 1, X-ray diffraction patterns obtained on blank 330Cb specimens oxidised for 48 h, in N 2 -5 vol.% H 2 at 900 °C indicate that after oxidation the scale is composed of Mn 1.5 Cr 1.5 O 4 (ICDD 33-0892) and Cr 2 O 3 (ICDD 38-1479). The intermetallic Nb 3 Ni 2 Si (ICDD 17-0097) is also detected after 48 h oxidation. This phase was not present on the alloy before oxidation and its presence is due to the silicon enrichment of the niobium nodules at high temperature. XRD patterns also show that, on the 7.10 -3 mg.cm -2 coated specimen the oxide scale is composed of Mn The scale cross-section morphology of the blank specimen oxidised in N 2 -5vol.% H 2 during 48 h show that a 1 µm thick oxide scale is present. Silicon is found associated with niobium inside the alloy and no silica layer is observed at the alloy/oxide interface. The scale cross-section morphology of the 7 10 -3 mg.cm -2 coated specimen oxidised in N 2 -5vol.% H 2 during 48 h is shown on Figure 2. The scale is mainly composed of a 2 µm thick manganese chromite scale. EDXS analysis shows that a continuous cristobalite scale is formed at the alloy/oxide interface. The sodium amount is too low to permit its detection in the scale. Oxides formed on the blank material, specimens coated with 0.007 mg.cm -2 , 0.068 mg.cm -2 and 0.262 mg.cm -2 NaOH. The scale cross-section morphology of the 6.8 10 -2 mg.cm -2 coated specimen oxidised in N 2 -5vol.% H 2 during 48 h is presented on Figure 3. After 48 hours oxidation at 900 °C in N 2 -5vol.% H 2 the cross-section micrograph shows that the scale is mainly composed of a 4 µm thick oxide scale. EDXS analysis shows that sodium and niobium are present together in the oxide scale indicating that sodium reacts with niobium provided by the alloy. A cristobalite layer is observed on the cross section.

The scale cross-section morphology of the 0.265 mg.cm -2 coated specimen oxidised in N 2 -5vol.% H 2 during 48 h is exhibited on Figure 4. The scale is mainly composed of a 6 µm thick iron chromite FeCr 2 O 4 scale. Chromium depletion is observed 10 µm deep inside the alloy. Sodium and niobium are present together in the oxide scale as a NaNbO 3 mixed oxide identified by XRD.

improves the oxidation resistance. Some authors stated that silicon segregates at the oxide/alloy interface and blocks the iron cationic diffusion [START_REF] Huntz | [END_REF][21][22]. Silicon is then supposed to be present as a silica film, which lowers the steel oxidation rate. The high silicon oxygen affinity permits its internal oxidation, developing SiO 2 precipitates at the internal interface [13,23]. Then, it is reported that silica acts as a diffusion barrier and also leads to the keying of the chromia scale to the substrate [24,25]. Silica will also lower the porosity at the internal interface acting as vacancies sinks [26]. Silicon reduces the amount of non-protective iron oxides inside the scale [27] and hinders the iron rich nodule formation [21]. It has been proposed that during the oxidation of a AISI 304 stainless steel between 900 et 1000 °C, a maximum 0.88 % (weight %) silicon content leads to the formation of a chromia scale, even at 1000 °C [28]. Nevertheless, it appears that high silicon content induces more scale spallation between the alloy and the silica scale or at the silica/chromia interface [22]. It is the reason why silicon is rarely added at more than 1 weight %. One should also take care that the chromium presence is necessary to avoid the fayalite Fe 2 SiO 4 formation, which acts as a very poor diffusion barrier [29]. According to Stott, the necessary amount of silicon needed to the SiO 2 formation can be lower when the chromium content increases [29]. Concerning the effect of higher silicon additions, Li and Gleeson compared the oxidation behaviour at 1000 °C of Ni-based alloys with and without about 2.7 wt. % Si additions [30]. From oxidation results of the cast model alloys, Si addition was found to improve oxidation resistance by forming a continuous SiO 2 layer at the alloy/scale interface, which resulted in decreased oxidation kinetics. The cast alloys, with silicon addition, also showed larger average effective inter-diffusion coefficients of chromium compared to the cast alloys without silicon addition. As a consequence, the silicon addition assisted in the establishment and re-formation of a chromia scale during oxidation. In the case of commercial alloys a discontinuous distribution of SiO 2 precipitates in the vicinity of the alloy/scale interface was found to be beneficial to cyclic oxidation resistance. The formation of an amorphous silica inner layer was effective in reducing the oxidation rate of a Fe-9%Cr alloy to about half that of Si-free steel (0.06%Si) at 700 °C, in steam, when alloyed with 0.5% Si. But the effect of silicon is less beneficial at lower temperatures (500 °C) [31]. Ishitsuda et al. carried out a detailed study on the effect of Si additions between 0.06% and 0.49% on the 9Cr steel in the temperature range 500 °C-700 °C in high pressure (35 MPa) steam, and concluded that silicon was most effective in reducing the oxidation rate at the highest temperature tested and at the highest silicon content, where rates were reduced to 50% of that observed with the low silicon alloy [32]. As with the work of Ueda et al. the establishment of a complete silica layer was responsible for the protective behaviour at high test temperatures [33,34].

In the present work, the alloy silicon content is 1.55 wt.%. Before oxidation, niobium nodules are present in the alloy and silicon is randomly distributed in the steel matrix. After oxidation of blank specimens EDXS results show that silicon is present in non-oxidised nodules at the scale/alloy interface and inside the alloy. Then, silicon is associated with niobium and nickel in the Nb 3 Ni 2 Si phase identified by XRD and no silica is formed.

The niobium content is 0.97 wt.% in the 330Cb alloy and should be taken into account. It is mentioned in the literature that, on austenitic steels (oxidised at 950 °C, 100h), niobium can favour the chromium diffusion in the metallic substrate to the steel/scale interface [35]. It has also been shown that the steel structures are more stable when niobium is associated with titanium and molybdenum [36]. Our results show the presence of Nb 3 Ni 2 Si nodules after oxidation of blank specimens. These nodules are sometimes close to the substrate surface and are not oxidised. The presence of such nodules can be explained by the fact that niobium only stays in solid solution in the alloy for temperature lower than 500 °C but not at higher temperatures [37,38]. Our results also show that when silicon is associated with niobium in Nb 3 Ni 2 Si intermetallic, niobium traps silicon during oxidation and hinders the silica scale formation.

On the other hand, XRD results indicate that silicon has been oxidised and cristobalite is present on sodium hydroxide coated specimens. SEM cross sections also indicate that silicon has diffused from the substrate to permit the formation of a silica scale located at the oxide/alloy interface (Figure 1 and Figure 2). The effect of a low oxygen environment, at the alloy surface, can be related to other works describing the effect of surface coatings limiting the oxygen access to the alloy surface [8]. Some authors have shown that an yttrium sol-gel coating promotes a continuous silica scale formation at the internal interface on the alloy AISI 304 due to a limited access of the oxygen to the metallic surface [39]. Then, silica hinders the iron oxidation and the formation of non-protective iron oxides [40]. Even though the alloy 330Cb contains 1.55 wt.% silicon, results show that no continuous silica scale formation was observed after oxidation in air, at 900 °C. Some scale spallation, observed on uncoated specimens, has also been explained by the presence of the high silicon amount in the alloy [41]. The absence of iron oxides in the scale can also be due to the fact that the alloy 330Cb is a nickel-rich alloy. Douglass [42] studied the Ni-2.05Si and Ni-4.45Si oxidation in oxygen for 18 hr in the temperature range 600-1000 °C. He proposed that the low silicon content alloy exhibited extensive internal oxidation and that the higher-silicon alloy formed a continuous layer of silicon-rich oxide. Nakakubo worked on Fe-Si and Fe-Cr alloys at 850 °C. This author calculated the boundary conditions between internal and external oxidation on the basis of the kinetic theory of the internal oxidation [43]. Onishi has also calculated the conditions at the boundary between internal to external oxidation of Si containing steels (Fe-Si alloys) at 850 °C [44]. This author has demonstrated that under low oxygen potential a very thin SiO 2 layer forms on the alloy surface because the diffusion rate of silicon is relatively high compared with that of oxygen. It is then concluded that internal oxidation occurs in a restricted oxygen partial pressure range and is limited to silicon concentrations below 1 wt.%. The oxygen partial pressure range, in which internal oxidation can occur, decreases with increasing silicon content. When the silicon content increases, external silicon oxidation should occur. In the present work, it is observed that under low oxygen potentials, the high silicon concentration in the 330Cb alloy permits the initial formation of a cristobalite scale on the alloys surface because the diffusion rate of silicon is relatively high compared with that of oxygen.

On 330Cb coated specimens, the coating was characterized by XRD and it was found that sodium hydroxide reacts with carbon dioxide to form sodium carbonate Na 2 CO 3 (ICDD 18-1208). Na 2 CO 3 decomposes at 856 °C to form Na 2 O and this compound is the one expected to react with the alloy at 900 °C. Results shown on figure 2 and figure 3 indicate that after oxidation of 7.10 -3 and 6.8 10 -2 mg.cm -2 coated specimens, a SiO 2 cristobalite scale is formed. XRD patterns show that when sodium hydroxide is deposited on the initial surface a large amount of NaNbO 3 phase could be detected in the oxide scale owing to its high amount (Figure 1). It should be remembered that when a crystalline phase is less present than 10 wt.%, it is hardly detected by XRD. SEM cross section (Figure 4) also shows that niobium and sodium are located at the same place. Due to the fact that the gaseous environment is always the same in this work, it is clear that sodium promotes the formation of cristobalite when low amounts of coatings are present (7.10 -3 and 6.8 10 -2 mg.cm -2 ). Results also indicate that no intermetallic Nb 3 Ni 2 Si could be detected by XRD after oxidation even though the alloy is still analysed. Moreover, XRD patterns shown on Figure 1 indicate that sodium reacts with niobium to form NaNbO 3 .

The NaNbO 3 formation leads to niobium depletion in the first 10 microns deep inside the alloy and not enough niobium is present to trap silicon in the Nb 3 Ni 2 Si intermetallic form. As a result, silicon is free to diffuse from the alloy matrix to the oxide/metal interface to form a silica scale as shown on figure 5. Results also indicate that low amounts of sodium hydroxide promote the silica crystallisation as proposed by other authors [45,46]. The cristobalite scale then acts as a good diffusion barrier against outward iron and chromium [40]. When a very well established cristobalite scale is built, a thin adherent oxide scale is formed and mainly composed of manganese chromite. The good oxide scale adherence is then related to the direct contact between SiO 2 cristobalite and manganese chromite and their similar dilatation coefficients, respectively 7 10 -6 K -1 and 7.5 10 -6 K -1 . When sodium hydroxide coatings exceed 0.265 mg.cm -2 , no cristobalite protective scale forms and iron combines with chromium to form FeCr 2 O 4 . As proposed by other authors high amounts of sodium salts can induce accelerated oxidation and hot corrosion processes [9][10][11]. It is then concluded that if low sodium additions are beneficial for the oxidation protection of the 330Cb alloy it should not exceed 6.8 10 -2 mg.cm -2 on the metallic surface.

CONCLUSION

This work shows the influence of sodium hydroxide coatings on the austenitic 330Cb (Fe-34Ni-23Cr-1Nb-1.55Si) oxidation during 48 h, at 900 °C. The N 2 -5vol.% H 2 gaseous environment was used in order to simulate industrial heat treatment conditions. With low NaOH deposits, the oxide scale is adherent and a SiO 2 cristobalite subscale is formed at the alloy/oxide interface. It is demonstrated that sodium combines with niobium to form NaNbO 3. As a niobium depleted area is formed in the alloy, results show that Nb 3 Ni 2 Si cannot be formed close to the internal interface and silicon is free to diffuse and oxidised as a silica scale at the oxide/alloy interface. On this austenitic steel, the protective silica scale formation is promoted by low oxygen partial pressures and the high silicon content. As the NaOH coating amount increases on the alloy surface, no protective silica scale is formed and severe oxidation is observed leading to fast growing FeCr 2 O 4 oxide scales. To avoid any catastrophic oxidation, the maximum amount of NaOH allowed on the alloy surface is about 6.8 10 -2 mg.cm -2 .

1 . 5

 15 Cr 1.5 O 4 (ICDD 33-0892) and Cr 2 O 3 (ICDD 38-1479). SiO 2 cristobalite (ICDD 39-1425) is also well detected on the specimen surface. The intermetallic Nb 3 Ni 2 Si (ICDD 17-0097) is not detected. On the 6.8.10 -2 mg.cm -2 coated specimen the oxide scale is composed of Mn 1.5 Cr 1.5 O 4 (ICDD 33-0892) and a higher amount of Cr 2 O 3 (ICDD 38-1479) is present on the surface. SiO 2 cristobalite (ICDD 39-1425) is present in lower concentration on the specimen surface compared to the previous one. On the 0.265 mg.cm -2 coated specimen, iron combines with chromium to form FeCr 2 O 4 (ICDD 34-0140). The proportion of Cr 2 O 3 (ICDD 38-1479) is low and sodium is found combined with niobium in a mixed oxide NaNbO 3 (ICDD 38-0846).

FIGURE 1 .

 1 FIGURE 1. XRD patterns obtained on the 330Cb alloy after 48 h oxidation at 900°C, in N 2 -5%H 2 .Oxides formed on the blank material, specimens coated with 0.007 mg.cm -2 , 0.068 mg.cm -2 and 0.262 mg.cm -2 NaOH.

FIGURE 2 .

 2 FIGURE 2. SEM cross sections obtained on the NaOH coated 330Cb specimen (7 10 -3 mg.cm -2 ) oxidized 48 h at 900 °C in N 2 -5vol.% H 2 (BSE Image x 5000).

3 FIGURE 5 .

 35 FIGURE 5. Schematic drawing showing the role of sodium-niobium interaction on the SiO 2 scale formation.

TABLE 1 . Alloy composition (weight %). Wt. %

 1

		Fe	Ni	Cr	Si	Nb	Mn	C	P	S
	330Cb	Bal.	34.41	22.87	1.55	0.97	0.72	0.05	0.012	0.002

Alloy

DISCUSSION