{'fraction_non_alphanumeric': 0.03618193465727098, 'fraction_numerical': 0.012581678411274824, 'mean_word_length': 4.428571428571429, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Accounting for missing actors in interaction network inference from abundance data

Network inference aims at unraveling the dependency structure relating jointly observed variables. Graphical models provide a general framework to distinguish between marginal and conditional dependency. Unobserved variables (missing actors) may induce apparent conditional dependencies.In the context of count data, we introduce a mixture of Poisson log-normal distributions with tree-shaped graphical models, to recover the dependency structure, including missing actors. We design a variational EM algorithm and assess its performance on synthetic data. We demonstrate the ability of our approach to recover environmental drivers on two ecological datasets. The corresponding R package is available from github.com/Rmomal/nestor.


Introduction
Network inference. Network inference (or structure inference) has become a topical problem in various fields such as biology, ecology, neuro-sciences, social sciences, to name a few. The aim is to unravel the dependency structure that relates a series of variables that can be jointly observed. Graphical models (see e.g. Lauritzen, 1996) provide a natural framework to achieve this task as it allows to encode the dependency structure into a graph, the nodes of which are the variables. Two variables are connected if and only if they are dependant, conditionally on all others. Most methodologies build on the assumption that the network is sparse, meaning that only a small fraction of variable pairs are conditionally dependent. The case of Gaussian graphical models (GGM) is especially appealing as the network corresponds to the support of the precision matrix of the joint Gaussian distribution. The use of a sparsity-inducing penalisation gives raise to the celebrated graphical lasso (Friedman et al., 2008). In a more general context, Chow and Liu (1968) consider a spanning tree structure to impose sparsity to the network, but this drastic form can be alleviated using mixtures of trees (Meilȃ and Jaakkola, 2006;Kirshner, 2008). One important aspect of network inference is to distinguish between variables that are marginally dependent (possibly because of their respective dependency with some common other) from variables that are directly related, that is conditionally dependant. This distinction requires to account for as many confounding effects as possible, which includes all the other variables but also available covariates. It also requires to consider the existence of some missing actors (or missing nodes), that may induce an apparent direct dependency.
Abundance data. Count data is found in a multitude of fields (sociology, biology, economy, ecology ...). It results from the counting of events in a given setting such as crime statistics in a state or the number of produced transcripts of a gene in an experiment. The statistical processing of count data cannot always rely on classical methods developed for continuous Gaussian data and appeals for specific methods. It often exhibits specificities such as zero-inflation and a large dispersion. The present work is motivated by the analysis of so-called abundance data, a count data avatar, arising from ecological studies where the number of individuals (the abundance) of a series of living species (plants or animals) is observed in a series of sites. In this context, network inference aims at understanding which pairs of species are in direct interaction. The covariates are typically environmental descriptors (altitude, temperature, distance to the see, etc.) of each collection site, while the variables are the respective abundances of each species from the community under study. No nice and generic framework as the GGM exists for count data. A few alternatives rely on copulas (Inouye et al., 2017) or models the node-wise conditional distributions as arising from exponential families. But most joint species distribution models resort to a latent Gaussian layer, which encodes the dependency structure between the species (Warton et al., 2015;Popovic et al., 2018Popovic et al., , 2019. The Poisson log-normal model (PLN: Aitchison and Ho, 1989) enters this category: it assumes that a multivariate Gaussian random variable is associated to each species in each site and that the observed abundances are conditionally independent Poisson variables. The PLN model has already been applied to abundance data, both for dimension reduction (Chiquet et al., 2018) and network inference (Chiquet et al., 2019;Momal et al., 2020).
Missing actors. In many situations, it is likely that not all actors involved in the system have been observed. The term 'actors' refers to either species that were not observed but nonetheless influence the abundance of others, or environmental conditions that were not accounted for.
In the perspective of unravelling the conditional independence structure, this can typically lead to the inference of spurious edges, which are links between observed actors that are not in direct interaction. In the graphical model framework, not accounting for one variable amounts to consider the marginal distribution of the rest of the system, as described in the left panel of Figure 1. Missing actors may be quantitative or qualitative. In the latter case it defines a latent group structure (Ambroise et al., 2009). Several approaches have been proposed for network inference accounting for quantitative missing actors in the context of GGM. Many of them (Chandrasekaran et al., 2011;Meng et al., 2014;Giraud and Tsybakov, 2012;Lauritzen and Meinshausen, 2012) adapted the principle of Robust PCA (Candès et al., 2011) to the concentration matrix, assuming it is a sum of two matrices: one low-rank and one sparse. In terms of missing actors in a network, the low-rank part corresponds to missing actors connected to all variables, whereas the sparse part refers to missing actors having a local effect. Following Robin et al. (2019) (also in the context of GGM), we focus on the later aspect, that is looking for missing actors not necessarily linked to all others. As far as we know, no model has been proposed for the inference of missing actors from abundance data.
Variational inference. The model we consider in this paper involves different types of variables, namely an unknown tree-shaped graphical model, a continuous latent layer (to induce dependence between the species) and unobserved actors. The most popular approach for the inference of such models is the EM algorithm (Dempster et al., 1977), which requires the evaluation of the conditional distribution of all unobserved variables given the data. In the problem we consider, some latent variables are (multivariate) continuous and others are discrete, and their joint conditional distribution turns out to be intractable. In this work we resort to a variational approximation (Wainwright and Jordan, 2008) of this conditional distribution and to a variational EM algorithm for its inference (see e.g. Blei et al., 2017).
Our contribution. In the context of the Poisson log-normal model, we propose a tree-based approach to recover the structure of latent graphical model including actors. The model we consider involves several layers of unobserved variables with intractable conditional distributions, thus we resort to a variational EM algorithm (Blei et al., 2017) for its inference. We introduce the model in Section 2 and describe its variational inference in Section 3. The performance of the algorithm is assessed via simulations in Section 4. The use of the proposed model is illustrated in Section 5, where we demonstrate its ability to recover environmental drivers on two ecological datasets. The inference procedure is implemented in the R package nestor, available at github.com/Rmomal/nestor.

Poisson log-normal and tree-shaped graphical models
Poisson log-normal model.
We start with a reminder on the multivariate Poisson log-normal model, with the example of abundance data. The abundances of p species observed on n sites are gathered in the n × p matrix Y where Y ij is the count of species j in site i, and the row i of Y , denoted Y i , is the abundance vector collected on site i. A covariate vector x i with dimension d is also measured on each site i and all covariates are gathered in the n × d matrix X. The PLN model states that a (latent) Gaussian vector U i of size p with variance matrix R = (ρ kl ) kl is associated to each site: the sites being assumed to be independent. To ensure identifiability, we let the diagonal of R be made of 1's, so R is actually a correlation matrix. All latent vectors U i are gathered in the n × p matrix U . The PLN model further assumes that species abundances in all sites are conditionally independent, and that their respective distribution only depends on the environment and the associated latent variable: where o ij is a known offset term which typically accounts for the sampling effort, and σ j is the latent standard deviation associated with species j. The vector d × 1 of regression coefficients θ j describes the environmental effects on species j. An important feature of the PLN model is that the sign of the correlation between the observed counts is the same as this of correlation between the latent variables (Aitchison and Ho, 1989): sign(Cor(Y ij , Y ik )) = sign(Cor(U ij , U ik )).

Tree-shaped graphical models.
Network inference relies on the assumption that few species are directly dependent on one another, meaning that the underlying graphical model is sparse. In the framework of the PLN model, the graphical model of interest rules the distribution of the latent vectors U i and is encoded in the precision matrix Ω := R −1 . A way to foster sparsity is to impose Ω to be faithful to a spanning tree T , that is: U 1 ∼ N p (0, Ω −1 T ) where the non-zero terms of Ω T correspond to the edges of the tree T . However this hypothesis is very restrictive as it allows only p − 1 links among p species (Chow and Liu, 1968). A more flexible approach consists in assuming that the latent vectors are drawn from a mixture of Gaussian distributions, each faithful to a tree T (Meilȃ and Jordan, 2000;Meilȃ and Jaakkola, 2006;Kirshner, 2008;Schwaller et al., 2019): where T p is the set of spanning trees with p nodes. We further assume that the tree distribution {p(T )} T ∈Tp can be written as a product over the edges: The weights β jk are gathered in the p × p symmetric matrix β with diagonal zero. Observe that these weights are defined up to a multiplicative constant, so that only p(p − 1)/2 − 1 of them may vary independently. This PLN model with latent tree-shaped dependency structure is similar to that considered by Momal et al. (2020).

Introducing the missing actor
PLN model with missing actors.
We now introduce the concept of missing actors, which corresponds to variables that are involved in the graphical model but are not associated to observed variables. To involve such actors in the model, we assume that a complete latent vector U i with dimension p + r is associated to site i, where r is the number of missing actors. This complete vector can be decomposed as where U Oi (with dimension p) corresponds to observed species and U Hi (with dimension r) corresponds to the missing actors. The complete n × (p + r) latent matrix U can be decomposed in the same way as U = [U O U H ], U O and U H having dimension n × p and n × r, respectively. The model we consider states that (i) the complete latent vectors U i are all iid and distributed according to a mixture similar to (3) and (4) but with Gaussian distributions (and matrices Ω T and β) of dimension (p + r), and trees drawn from T p+r ; (ii) the abundances Y ij of the p observed species are distributed according to (2) In the sequel, we shall refer to the elements of U O and U H respectively as 'observed' and 'hidden' (or 'missing') latent variables, whereas obviously none of them are actually observed. Figure 2 displays the graphical model of the quadruplet (T, U O , U H , Y ). The observed data Y still arise from an PLN model, but the graphical model of the observed latent U O may not be sparse due to the marginalization over the hidden latent U H . Our main goal is to infer the dependency structure of the complete latent vectors, that is to estimate the elements of the matrices Ω T and the edges weights β. The latent dependency structure is similar to this considered by Robin et al. (2019), but the inference strategy much differs, because of the additional hidden layer.

Identifiability restriction.
The proposed model only makes sense because the graphical model of the complete latent vectors U i = [U Oi U Hi ] is supposed to be sparse. Missing actors could obviously not be identified from a regular PLN model, without restriction on the precision matrix Ω, as only the marginal precision matrix of the U Oi could be recovered. Still, to ensure identifiability we impose the same restriction as Robin et al. (2019) that missing latent variables are not connected with each other (the block corresponding to U H × U H is diagonal in each Ω T ).

Inference
As said in the introduction, we resort to a variational EM algorithm to perform the inference due to the complex latent structure.

Variational inference
The log-likelihood of the so-called complete data, that is (Y , U , T ), writes where Ω stands for the set of all tree-specific precision matrices: Ω = {Ω T , T ∈ T p+r }. The conditional distributions of the latent variables U and of the tree T given the data Y are both intractable. Variational inference then aims at maximizing a lower bound of the log-likelihood of the observed data, which writes in our context as where q(U , T ) stands for the approximate joint conditional distribution of the latent layer and of the tree: q(U , T ) p(U , T | Y ).
Approximate distribution.
The efficiency of variational inference mostly depends on the choice of q(U , T ), which is a balance between computational ease and adequation to the target distribution p(U , T | Y ). We adopt here a classical product form for the approximate distribution: we impose to the latent variables U and to the tree T to be independent according to q (whereas actually they are not conditional on the data), with respective marginals h and g: Because the sites are independent, and without further assumption, the distribution h is a product over all sites. Following Chiquet et al. (2018) we approximate the conditional distribution of each latent vector U i with a Gaussian distribution, that is: with all S i diagonal. We gather all the mean vectors m i in the n × (p + r) matrix M and pile up the diagonals of all the variance matrices S i in the n × (p + r) matrix denoted S.

Variational EM.
The variational EM algorithm then consists in maximizing the lower bound J defined in (5) with respect to the parameters (M step), and to the approximate distributions (VE step), alternatively.
M step: At iteration t + 1, given the current approximate distribution q t (U , T ) = g t (T )h t (U ), the M step consists in the update of the model parameters, solving Observe that the matrix of edge weights β is considered here as a parameter to be estimated, as opposed to Robin et al. (2019), where is was kept fixed and supposed to be given.

VE step:
Maximising J with respect to (wrt) q is equivalent to minimizing the Küllback-Leibler divergence between q(U , T ) and p θ,β,Ω (U , T | Y ) that appears in (5). Because we adopted a product form for q, the solution of the VE step for both g and h is known to be a mean-field approximation (Wainwright and Jordan, 2008). More specifically, maximising J gives and Observing that log p β (T ) + log p Ω (U | T ) can be written as a sum over all the edges present in T , we see that g t+1 (T ) has a product form. So, without any further assumption, we may parametrize g(T ) in the same way as p β (T ): We gather the β jk 's in the (p + r) × (p + r) matrix β. The parameters β, M and S are called the variational parameters, in the sense that it is equivalent to optimize J wrt (g, h) or wrt ( β, M , S).

Proposed algorithm
The model we consider is an extension of the PLN model, for which an efficient inference algorithm have been implemented in the PLNmodels, an R package available on CRAN (Chiquet et al., 2018(Chiquet et al., , 2019.

Prior estimates of θ, M O and S O .
To alleviate the computational burden of the inference, we take advantage of this available tool to get an estimate of the regression coefficient matrix θ and an approximation of the parameters of the observed latent variable conditional distribution  This steps deals with the update of the model parameters β and Ω T . Some of the calculations are tedious and postponed to Appendix B. (6), the maximization of J requires the computation of the derivative of E g t [log p β (T )] wrt β, which includes the derivative of the normalizing constant B. The latter can be computed via an extension of the Matrix Tree theorem (see Meilȃ and Jaakkola, 2006, Lemma 1 reminded in Appendix A). Setting the derivative of the expectation to 0 yields the following update (same as in Momal et al. (2020) and detailed in appendix B.1):

Edges weights β: As shown in Equation
where M (β) is defined in Lemma 1 and P t kl is the probability that the edge (k, l) belongs to the tree T according to g t : P t kl is computed using a result from Kirshner (2008) (reminded as Lemma 2 in appendix A). We now define the binary variable I T kl which indicates the presence of the edge kl in tree T , so P t kl = E g t [I T kl ] and I T = [I T kl ] 1≤k,l≤(p+r) is the adjacency matrix of tree T .
Precision matrices Ω T : For a given dependency structure in the Gaussian Graphical model framework, Lauritzen (1996) gives maximum likelihood estimates for the precision matrix. These estimators are given as functions of sufficient statistics of the multivariate Gaussian distribution. Indeed in the exponential family framework, the M step of any EM algorithm requires the computation of the expectation of a sufficient statistic, under the current fit of the variational laws (see McLachlan and Krishnan (2007)). Here as U | T is centered, a sufficient statistic is U U . We now let SSD denote the matrix defined as Applying Lauritzen's formulas, we get: where ssd t kl stands for the entry kl of the matrix SSD t . The calculations are postponed to Appendix B.2. Observe that estimates of the off-diagonal entries ω t+1 T kl do not depend on T provided that the edge (k, l) belongs to T . Thus the estimates of the off-diagonal terms of the precision matrices Ω T are common to all trees sharing a given edge. This does not result from any assumption on the shape of Ω T , but from the properties of the maximum likelihood estimate of Gaussian variance matrix. In the sequel we will simply denote off-diagonal terms by ω kl (as opposed to ω T kk which still depends on T ).
Other quantities are needed for later computations. Lauritzen gives the maximum likelihood estimator of every entry of the correlation matrix R T corresponding to an edge kl being part of T , which is R t+1 T kl = ssd t kl /n. Hereafter for any matrix A, A [kl] refers to the bloc kl of A: A [kl] = (a ij ) {i,j}∈{k,l} . The determinant of Ω t+1 T factorizes on the edges of T and writes as a function of blocs of the correlation matrix as follows: Finally we define the matrix Ω kl , edges probabilities appear as follows:

VE step.
This step deals with the update of the approximate conditional distributions g and h H , namely the update of the corresponding variational parameters β, M H and S H .
Approximate conditional tree distribution g(T ): Computing the expression (7) yields the following, where the constant term 'cst' does not depend on a specific edge: Then remembering the product form of g t+1 given in (9), we obtain the expression for each edge variational weight: Approximate Gaussian distribution h: According to (8), we have that Using the properties of the conditional Gaussian distribution we have that , Ω
As for any EM algorithm, the choice of the starting point is paramount. The initialization we use here takes the primary estimate M O as an input.
Initial clique: As a starting point, we look for a clique of species as potential neighbors of the missing actor h. There are many different ways to do so, and if any prior knowledge exists on that matter it should be used. Otherwise, such a clique can be found using sparse principal component analysis (sPCA; Erichson et al., 2020), where principal components are formed using only a few of the original variables, which is consistent with the assumption that each missing actor is connected only to some actors in the network.
When applying sPCA to M O , the set of non-zero loadings of each principal components provides us with an initial clique of neighbors of each missing actor Parameters initialization: The eigenvectors resulting from the sPCA also provide us with a starting value M 0 H , as well as a first estimate of the latent correlation matrix R 0 . The parameter β is uniformly initialized.

Numerical issues.
Because the Matrix Tree Theorem and Kirshner's formula respectively resort to the calculation of a determinant and a matrix inversion, the proposed algorithm is exposed to numerical instabilities. To circumvent these issues, we rely on both multiple-precision arithmetic and likelihood tempering (via a parameter α, similarly to Schwaller and Robin, 2017). More details are given in Appendix B.4.

Count datasets
For the simulation study, 300 count datasets of 15 species in total including one missing actor are generated, thus p = 14 and r = 1. Data is generated as follows. We generate a scale-free structure G (which degree distribution is a power law) with p + 1 nodes using the R package huge (Zhao et al., 2012) available on CRAN. The missing species h is chosen as the one with highest degree. We measure the influence of the missing actor with its degree, distinguishing three influence classes: Minor (degree ≤ 5), Medium (5 < degree ≤ 7) and Major (degree ≥ 8). For each replicate, the latent layer U and the observed abundances Y are simulated according to the model defined in Section 2.

Experiment & Measures
For each simulated dataset, the VEM algorithm is initialized as described in Section 3.3. More specifically and because we only look for one missing actor, we consider the cliques corresponding to each of the first two principal components of sPCA, and their respective complements, which provides us with four cliques. Then four VEM algorithms, as described in Section 3.2, are run starting from each of the four candidate cliques, and the one yielding the highest lower bound J is kept. For all simulations, we set the precision of the convergence criterion to ε = 10 −3 , the tempering parameter to α = 0.1 and the maximal number of iterations to 100. The inference quality is assessed regarding the global network inference, the missing actor's position in the network, and its values along the n sites. We refer to this first procedure as the blind procedure. Additionally, we define the oracle procedure as running the VEM with the set of true neighbors of the missing actor as initial clique.
For each procedure, a general measure of the whole network inference quality is first given by comparing the inferred edge probabilities to the original dependency structure. This is done using the Area Under the ROC Curve (AUC) criteria. Then, to be more specific and target the neighbors of node h specifically, the probabilities of edges involving h are transformed into binary values using the 0.5 threshold. The values are then compared to the original links of h and yield quantities of true/false (T /P ) positives/negatives (P /N ), from which are built the precision (also known as the positive predictive value, T P /(T P + F P )) and the recall (also known as the true positive rate, T P /(T P + F N )) criteria. Finally, we assess the ability to reconstruct the missing actor across the sites by computing the absolute correlation between its inferred vector of means (M h ) and its original latent Gaussian vector U h .

Results
Simulations performance measures are gathered in Table 1 and Table 2 for blind and oracle procedures respectively. The distributions of the quality measures are displayed in Figure 3. Table 1 shows the network is well inferred, as all AUC means are above 0.85, with almost perfect inference when the influence of the missing actor is major. Its neighbors and values per site are very well retrieved in these cases with mean recall values above 0.9 and mean correlation above 0.8, with a great confidence in the algorithm outputs as mean precision is above 0.95. However, there exists a clear deterioration of all performance as the influence decreases with lower means are greater deviations, down to about 0.6 mean values for all measures when the influence is minor. Moreover, the algorithm takes more and more time to converge as the influence decreases, although it stays at about 3s for minor cases which is reasonable. Figure 3 shows that as the influence decreases, the densities present with several modes and dilute towards 0, illustrating that even if some networks are still well-inferred, there also are more and more cases where the algorithm fails. In particular, the performance decrease of medium cases seems to be only due to a greater number of failed inferences.
All these elements point to minor cases being harder problems to solve, unsurprisingly. Yet as oracle results show in Table 2, it is possible to carry out almost-perfect inference in all cases, if the algorithm is initialized with the true clique; the deterioration is still present in all measures, but stays marginal. Thus the harsh decrease in the blind procedures seems to be mainly due to the proposed initialization method failing at correctly finding some of the small cliques of neighbors.
About intialization. Figure 4 compares the initialization quality and the corresponding final inferred neighbors, in terms of initial (-i) and final (-f) false negative (FNR, also 1-TPR) and positive rates (FPR). It clearly appears that final measures mostly increase with false negatives of the initial clique. This means that not including a neighbor in the initialization is much worse for the inference than falsely including a node. The increase of FNR-f is bigger than that of FPR-f, meaning that a wrong initialization leads to a set of inferred neighbors which most part can be trusted, but which will be largely incomplete. This advocates for bigger initialization cliques when no prior information is available.  Table 2: Oracle procedure using true clique as starting point. The influence of the missing actor is measured with its degree, distinguishing three influence classes: Minor (degree ≤ 5), Medium (5 < degree ≤ 7) and Major (degree ≥ 8). For each class of influence, the following quantities are reported: number of simulated graphs (N), means and standard deviations of AUC, Precision, Recall, Correlation between missing actor inferred vector of means and original latent vector, and running times in seconds. AUC measures the retrieval of the dependence structure between all variables (observed and missing), whereas precision and recall are specific to the missing actor links.

Cross validation criterion for model selection
The proposed model obviously raises the problem of choosing the number of missing actors r (which may be zero). Variational-based inference often relies on approximate versions of the BIC or ICL criteria for model selection. Few theoretical guaranties exist about these approximate criteria and, in the present case, we observed that BIC and ICL penalizations did not yield consistent results. Therefore, we resort to V -fold cross validation to determine the number of missing actors. More specifically, we split the original dataset Y (X is dropped here for the sake of clarity) into V subsets with almost equal sizes m 1 , .
For each subset v, we define its complement Y −v on which we fit a model with r missing actors and get a parameter estimate To avoid the integration over the (p + r)-dimensional Gaussian latent layer, we measure the fit with the pairwise composite likelihood (Lindsay, 1988). For any given tree T and parameter Γ, the bivariate Poisson log-normal pdf p P LN ((Y ij , Y ik ); Γ, T ) can be easily computed for any sample i and pair of species (j, k) with available tools such as the poilog R package (Vidar and Steinar, 2008) available on CRAN. The cross-validation criterion is defined as where the tree samples {T −v r,b } b=1...B are iid according to p β −v r (T ). The sampling procedure for spanning trees is given in Appendix C.1; the complete procedure for the calculation of P CL r (Y ) is described by Algorithm 1, given in Appendix C. Note that this criterion measures the fit of the model in terms of abundance prediction, whereas our interest is mostly focused on the inference of the dependency structure. In other words, our goal is identification, that is selecting the smallest model and not the best model in terms of prediction (Arlot and Celisse, 2010).
We did not include this computationally greedy procedure in the simulation study but applied it to the two ecological datasets that will be described in the next two sections. The results, gathered in Figure 5, yield r = 1 missing actor for the Barents Sea data set, and r = 2 missing actors for the Fatala River one. Regarding the initialization, we performed a wider exploration as compared to the simulation study. To enlarge the list of possible cliques, we applied a resampling version of the procedure described in Section 3.3, and applied it to 200 sub-samples, each consisting in 80% of the whole data set. This yielded 200 lists of r initial cliques, from which duplicates were removed.

Barents Sea
The dataset was first published by Fossheim et al. (2006) and consists of the abundance of 30 fish species measured in 89 sites in the Barents See in April-May 1997. In addition to abundances, the water temperature was measured in each site. The complete dataset is available at www.fbbva.es/microsite/multivariate-statistics/data.html. Fishes distributions are known to be greatly linked with the temperature. Hence to illustrate our methodology, we present the results of the model fitted without any covariate (that is not accounting for the temperature), but including one missing actor (as suggested by Figure 5). To assess the ability of the proposed methodology to retrieve the influence of temperature as a missing actor, we report the empirical correlation between the temperature and the conditional expectation of the missing actor M h , which we denote ρ(H, temp).
The resampling initialization procedure yielded in 14 different cliques, for each of which a VEM algorithm was run: the mean running time was 6.63mins with deviation 0.70 mins.
The edge probabilities involving node h as an endpoint were either very close to 0 or very close to 1, yielding a total of 6 highly probable neighbors of h. Figure 6 shows that many direct interactions are inferred between the corresponding 6 species in absence of a missing actor, which vanish when it is introduced. It also shows that accounting for this actor has only a local effect and that the direct interactions among the other species are preserved, which is consistent with our notion of a missing actor.  In terms of interpretation, Figure 7 shows that the missing actor is highly correlated with the temperature. It also appears that the abundances of the species neighbor to the missing actor are much more correlated with the temperature (mean correlation = 0.78, sd = .06) than the abundances of the non-neighbor species (mean correlation = 0.46, sd = .27). This example shows the ability of the method to recover an underlying effect that would not be recorded in the data. Baran (1995) collected the abundances of 33 fish species in 90 sites along the Fatala River in Guinea between June 1993 and February 1994. The data are available from the R package ade4 on CRAN (Dray et al., 2007), along with the date and site of collection, from which we deduce the season (dry or rainy). Again the model was fitted without any covariates, but with two missing actors, as suggested by Figure 5. The resampling initialization procedure yielded in 60 different cliques, for each of which a VEM algorithm was run: the mean running time was 11.33 min (sd = 1.47 mn). 14 VEM did not reach convergence (with tolerance ε = 1e − 3) after 100 iterations. We filtered out the results obtained from the different initializations, when the algorithm obviously ended in a degenerate solution (V(M h ) < exp(−20)). Figure 8 shows the scatterplot of the estimated conditional mean of the two missing actors (M h1 , M h2 ) in each site, colored with either one of the available covariates (site and season). The missing actor h 1 is obviously linked to the site and separates most upstream locations (kilometer 3) from most downstream locations (kilometer 46). This actor has 11 highly probable neighbor species. Again, this retrieved missing actor corresponds to an underlying effect (in this case: geography) that rules fish species abundances. The second missing actor seems to be linked with the season but with a less clear separation. Also the variability of M h2 is much smaller than this of M h1 . This effect is therefore questionable, which brings us back to model selection. As mentioned above, we used a procedure based on cross-validation, which may be prone to select too complex model (Shao, 1993;Friedman et al., 2001;Arlot and Celisse, 2010). The definition of a grounded model selection criterion for structure inference in presence of missing actors remains open.

A Algebraic Tools
We here present some algebraic results about spanning tree structures which are used during the computations. Theorem 1, Lemma 1 as well as Lemma 2 use the notion of Laplacian matrix Q of a symmetric matrix W = [w jk ] 1≤j,k≤p , which is defined as follows : We further denote W uv the matrix W deprived from its uth row and vth column and we remind that the (u, v)-minor of W is the determinant of this deprived matrix, that is |W uv |.
The following Theorem 1 is the extension of Kirchhoff's Theorem to the case of weighted graphs (Chaiken and Kleitman, 1978;Meilȃ and Jaakkola, 2006).
Theorem 1 (Matrix Tree Theorem) For any symmetric weight matrix W with all positive entries, the sum over all spanning trees of the product of the weights of their edges is equal to any minor of its Laplacian. That is, for any 1 ≤ u, v ≤ p, In the following, without loss of generality, we will choose Q 11 . As an extension of this result, Meilȃ and Jaakkola (2006) provide a close form expression for the derivative of W with respect to each entry of W .
Lemma 1 (Meilȃ and Jaakkola (2006)) Define the entries of the symmetric matrix M as it then holds that Kirshner (2008) build on Lemma 1 to provide an efficient computation of all edges probabilities.
Lemma 2 (Kirshner (2008)) Let p W be a distribution on the space of spanning trees, such that p W (T ) = kl∈T w kl /W , where W is defined as in Theorem 1. Taking the symmetric matrix M as defined in Lemma 1, the probability for an edge kl to be in the tree T * writes: As in Momal et al. (2020), the update of β is such that: By definition of p β (T ): Computing the derivative with respect to the edge weight β kl gives: According to Lemma 1: ∂ β kl B t = [M ] kl ×B. Finally setting the derivative to 0 yields the update

B.2 Update of Ω T
The update of Ω T respects This is a problem of parameter optimisation in the context of Gaussian Graphical Models (GGM).
In what follows, for any q × q matrix A, A [kl] will refer to the bloc kl of A: A [kl] = (a ij ) {i,j}∈{k,l} .
[A [kl] ] q will then denote the matrix obtained by filling up with zero entries to obtain full dimension q × q, so that: In its proposition 5.9, Lauritzen (1996) states that in a GGM with p variables and associated with the decomposable graph G, the maximum likelihood of the precision matrix exists if and only if n > max C∈C |C|. It is then given as where C is the set of cliques and S the set of separators of G, with associated multiplicities ν(S).
In our context, G is a spanning tree and so all cliques are edges and separators are nodes. The multiplicity of a given node k as a separator in the graph is ν(k) = d(k) − 1, where d(k) is its degree. Therefore the estimator Ω T writes as the following As SSD has diagonal n, the expression simplifies. Denoting I d the identity matrix of dimension d we obtain: Ω T = n kl∈T [(SSD [kl] ) −1 − 1 n I 2 ] p+r + I p+r .
The determinant of a precision matrix of a GGM with a decomposable graph is expressed as follows (Lauritzen, 1996): As Ω T is tree-structured, its determinant factorizes on the edges of T . It is expressed with the correlation matrix R T as follows, denoting d(k) the degree of node k: Using that R T has diagonal 1, we obtain for step t + 1 of the algorithm:

Exact computations
Our algorithm requires the computation of determinants (from the Matrix Tree Theorem) and inverses (in Kirshner's formula) of Laplacian of weight matrices. As we deal with highly variable weights, numerical issues arise: infinite determinants or matrix numerically non-invertible due to either the maximal machine precision (about 1.7·10 308 ), or with machine zero (about 2.2·10 −16 ).
To enhance the precision of such computations, we rely on multiple-precision arithmetic which allows the digit of precision of numbers to be limited only by the available memory instead of 64 bits. We implemented matrix inversion and log-determinant computation using both, symbolic computation and multiple precision arithmetic, relying on the gmp R package available on CRAN, which uses (Lucas et al., 2020), the C library GMP (GNU Multiple Precision Arithmetic).

Tempering parameter α
Weights β are mechanically linked to the quantity of data available n. To avoid reaching maximal precision when computing the determinant, a tempering parameter α is applied to every quantity proportional to n, so that the actual update performed is log β kl = log β kl − α( n 2 log | R T kl | + ω T kl [M M ] kl ).
We provide hereafter a heuristic to set the parameter α.
An upper bound for α: The proposed algorithm requires the computation of the normalizing constant B, which is the determinant of any minor of the Laplacian of the q×q variational weights matrix β. As these weights mechanically increase with the quantity of available data n, this step is numerically very sensitive. Hereafter we denote |Q uv | this determinant and ∆ the maximal machine precision. In order to ease the computations, we define the tempering parameter α as log β kl = log β kl − α( n 2 log | R T kl | + ω T kl [M M ] kl ) , under constraint |Q uv | ≤ ∆.