Figure 1 shows a visual schematic of the FEEMS method. The input data are genotypes and spatial locations (e.g. latitudes and longitudes) for a set of individuals sampled across a geographic region. We construct a dense spatial grid embedded in geographic space where nodes represent sub-populations, and we assign individuals to nodes based on spatial proximity (see Appendix 1—figure 1 for a visualization of the grid construction and node assignment procedure). The density of the grid is user defined and must be explored to appropriately balance model mis-specification and computational burden. As the density of the lattice increases, the model is similar to discrete approximations used for continuous spatial processes, but the increased density comes at the cost of computational complexity.

Figure 1

Download asset Open asset

Details on the FEEMS model are described in the Materials and methods section, however at a high level, we assume exchangeability of individuals within each sub-population and estimate allele frequencies, ${\widehat{f}}_j\left(k\right)$, for each sub-population, indexed by $k$, and single nucleotide polymorphism (SNP), indexed by $j$, under a simple Binomial sampling model. We also use the recorded sample sizes at each node to model the precision of the estimated allele frequency. With the estimated allele frequencies in hand, we model the data at each SNP using an approximate Gaussian model whose covariance is, up to constant factors, shared across all SNPs—in other words, after rescaling by SNP-specific variation factors, we assume that the set of observed frequencies at each SNP is an independent realization of the same spatial process. The latent frequency variables, $f_j\left(k\right)$, are modeled as a Gaussian Markov Random Field (GMRF) with a sparse precision matrix determined by the graph Laplacian and a set of residual variances that vary across SNPs. The pseudo-inverse of the graph Laplacian in a GMRF is inherently connected to the notion of resistance distance in an electrical circuit (Hanks and Hooten, 2013) that is often used in population genetics to model the genetic differentiation between sub-populations (McRae, 2006). The graph’s weighted edges, denoted by $w_{ij}$ between nodes $i$ and $j$, represent gene-flow between the sub-populations (Friedman et al., 2008; Hanks and Hooten, 2013; Petkova et al., 2016). The Gaussian approximation has the advantage that we can analytically marginalize out the latent frequency variables. The resulting likelihood of the observed frequencies shares a number of similarities to that of EEMS (see Materials and methods).

To prevent over-fitting we use penalized maximum likelihood to estimate the edge weights of the graph. Our overall goal is thus to solve the following optimization problem:

where $\mathit{𝒘}$ is a vector that stores all the unique elements of the weighted adjacency matrix, $\mathit{𝒍}$ and $\mathit{𝒖}$ are element-wise non-negative lower and upper bounds for $\mathit{𝒘}$, $\mathrm{\ell }\left(\mathit{𝒘}\right)$ is the negative log-likelihood function that comes from the GMRF model described above, and ${\varphi }_{\lambda }\left(\mathit{𝒘}\right)$ is a penalty that controls how constant or smooth the output migration surface will be and is controlled by the hyperparameter $\lambda >0$. Writing $𝒱$ to denote the set of nodes in the graph and $ℰ\left(i\right)\subset 𝒱$ to denote the subset of nodes that have edges connected to node $i$, our penalty is given by

This function serves to penalize large differences between the weights $w_{ik}$ and $w_{i\mathrm{\ell }}$ on edges that are adjacent, that is, penalizing differences for any pair of edges that share a common node. The tuning parameter λ controls the overall strength of the penalization placed on the output of the migration surface—if λ is large, the fitted surface will favor a homogeneous set of inferred migration weights on the graph, while if λ is low, more flexible graphs can be fitted to recover richer local structure, but this suffers from the potential for over-fitting. The tuning parameter λ is selected by evaluating the model’s performance at predicting allele frequencies at held out locations using leave-one-out cross-validation (see Materials and methods ‘Leave-one-out cross-validation to select tuning parameters’).

The scale parameter ${\widehat{w}}_0$ is chosen first fitting a ‘constant $w$’ model, which is a spatially homogeneous isolation-by-distance model constrained to have a single $w$ value for all edges. In the ${\varphi }_{\lambda }$ penalty, for adjacent edges $(i,k)$ and $(i,\mathrm{\ell })$, if $w_{ik}$ and $w_{i\mathrm{\ell }}$ are large (relative to ${\widehat{w}}_0$) then the corresponding term of the penalty is approximately proportional to ${\left(w_{ik}-w_{i\mathrm{\ell }}\right)}^2$, penalizing differences among neighboring edges on a linear scale; if instead $w_{ik}$ and $w_{i\mathrm{\ell }}$ are small relative to ${\widehat{w}}_0$, then the penalty is approximately proportional to ${\left(\log \left(w_{ik}\right)-\log \left(w_{i\mathrm{\ell }}\right)\right)}^2$, penalizing differences on a logarithmic scale. In fact, it is also possible to consider treating this scale parameter as a second tuning parameter—we can define a penalty function ${\varphi }_{\lambda ,\alpha }\left(\mathit{𝒘}\right)=\frac{\lambda }{2}{\sum }_{i\in 𝒱}{\sum }_{k,\mathrm{\ell }\in ℰ\left(i\right)}{\left(\log \left(e^{\alpha w_{ik}}-1\right)-\log \left(e^{\alpha w_{i\mathrm{\ell }}}-1\right)\right)}^2$, and explore the solution across different values of both λ and α. However, we find that empirically choosing $\alpha =1/{\widehat{w}}_0$ offers good performance as well as an intuitive interpretation (i.e. scaling edge weights $w_{ik}$ with reference to the constant-$w$ model), and allows us to avoid the computational burden of searching a two-dimensional tuning parameter space.

We use sparse linear algebra routines to efficiently compute the objective function and gradient of our parameters, allowing for the use of widely applied quasi-Newton optimization algorithms (Nocedal and Wright, 2006) implemented in standard numerical computing libraries like scipy (Virtanen et al., 2020) (RRID:SCR_008058). See the Materials and methods section for a detailed description of the statistical models and algorithms used.

While our statistical model is not directly based on a population genetic process, it is useful to see how it performs on simulated data under the coalescent stepping stone model (Figure 2, also see Appendix 1—figure 2 for additional scenarios). In these simulations we know, by construction, the model we fit (FEEMS) is different from the true model we simulate data under (the coalescent), allowing us to assess the robustness of the fit to a controlled form of model mis-specification.

Figure 2

Download asset Open asset

The first migration scenario (Figure 2A–C) is a spatially homogeneous model where all the migration rates are set to be a constant value on the graph, this is equivalent to simulating data under an homogeneous isolation-by-distance model. In the second migration scenario (Figure 2D–E), we simulate a non-homogeneous process by representing a geographic barrier to migration, lowering the migration rates by a factor of 10 in the center of the habitat relative to the left and right regions of the graph. Finally, in the third migration scenario (Figure 2G–I), we simulate a pattern which corresponds to anisotropic migration with edges that point east/west being assigned to a fivefold higher migration rate than edges pointing north/south. For each migration scenario, we simulate two sampling designs. In the first ‘dense-sampling’ design (Figure 2B,E,I) we sample individuals for every node of the graph. Next, in the ‘sparse-sampling’ design (Figure 2C,F,J) we sample individuals for only a randomly selected 20% of the nodes.

For each coalescent simulation, we used leave-one-out cross-validation (at the level of sampled nodes) to select the smoothness parameter λ. In the homogeneous migration simulations, the best value for the smoothness parameter, as determined by the grid value with the lowest leave-one-out cross-validation error, is ${\lambda }_{\text{cv}}=100$ in both sampling scenarios with complete and missing data. In the heterogeneous migration simulations ${\lambda }_{\text{cv}}=0.298$ with no missing data and ${\lambda }_{\text{cv}}=37.927$ with missing data. Finally, in the anisotropic simulations with no missing data ${\lambda }_{\text{cv}}=0.298$ and with missing data ${\lambda }_{\text{cv}}=0.042$. We note the magnitude of the selected λ depends on the scale of the loss function so comparisons across different datasets are not generally interpretable.

With regard to the visualizations of effective migration, FEEMS performs best when all the nodes are sampled on the graph, that is, when there is no missing data (Figure 2B,E,H). Interestingly, in the simulated scenarios with many missing nodes, FEEMS can still partly recover the migration history, including the presence of anisotropic migration (Figure 2I). A sampling scheme with a central gap leads to a slightly narrower barrier in the heterogeneous migration scenario (Appendix 1—figure 2I) and for the anisotropic scenario, a degree of over-smoothness in the northern and southern regions of the center of the graph (Appendix 1—figure 2N). For the missing at random sampling design, FEEMS is able to recover the relative edge weights surprisingly well for all scenarios, with the inference being the most challenging when there is anisotropic migration. The potential for FEEMS to recover anisotropic migration is novel relative to EEMS, which was parameterized for fitting non-stationary isotropic migration histories and produces banding patterns perpendicular to the axis of migration when applied to data from anisotropic coalescent simulations (Petkova et al., 2016, Supplementary Figure 2; see also Appendix 1 ‘Edge versus node parameterization’ for a related discussion). Overall, even with sparsely sampled graphs, FEEMS is able to produce visualizations that qualitatively capture the migration history in coalescent simulations.

To assess the performance of FEEMS on real data, we used a previously published dataset of 111 gray wolves sampled across North America typed at 17,729 SNPs (Schweizer et al., 2016; Appendix 1—figure 5). This dataset has a number of advantageous features that make it a useful test case for evaluating FEEMS: (1) The broad sampling range across North America includes a number of relevant geographic features that, a priori, could conceivably lead to restricted gene-flow averaged throughout the population history. These geographic features include mountain ranges, lakes, and islands. (2) The scale of the data is consistent with many studies for non-model systems whose spatial population structure is of interest. For instance, the relatively sparse sampling leads to a challenging statistical problem where there is the potential for many unobserved nodes (sub-populations), depending the density of the grid chosen.

Before applying FEEMS, we confirmed a signature of spatial structure in the data through regressing genetic distances on geographic distances and top genetic PCs against geographic coordinates (Appendix 1—figures 6, 7, 8, 9). We also ran multiple replicates of ADMIXTURE for $K=2$ to $K=8$, selecting for each $K$ the highest likelihood run among replicates to visualize (Appendix 1—figure 10). As expected in a spatial genetic dataset, nearby samples have similar admixture proportions and continuous gradients of changing ancestries are spread throughout the map (Bradburd et al., 2018). Whether such gradients in admixture coefficients are due to isolation by distance or specific geographic features that enhance or diminish the levels of genetic differentiation is an interpretive challenge. Explicitly modeling the spatial locations and genetic distance jointly using a method like EEMS or FEEMS is exactly designed to explore these types of questions in the data (Petkova, 2013; Petkova et al., 2016).

We first show FEEMS results for four different values of the smoothness parameter, λ from large $\lambda =100$ to small ${\displaystyle \lambda =0.0008}$ (Figure 3). One interpretation of our regularization penalty is that it encourages fitting models of homogeneous and isotropic migration. When λ is very large (Figure 3A), we see FEEMS fits a model where all of the edge weights on the graph nearly equal the mean value, hence all the edge weights are colored white in the relative log-scale. In this case, FEEMS is fitting a relatively homogeneous migration model where all the estimated edge weights get assigned nearly the same value on the graph. As we sequentially lower the penalty parameter, (Figure 3B,C,D) the fitted graph begins to appear more complex and heterogeneous as expected (discussed further below). Figure 3E shows the cross-validation error for a pre-defined grid of λ values (also see Appendix 1—figure 6 for visualizations of the fitted versus genetic distance on the full dataset).

Figure 3

Download asset Open asset

The cross-validation approach finds the optimal value of λ to be 2.06. This solution visually appears to have a moderate level of regularization and aligns with several known landscape features (Figure 4). Spatial features in the FEEMS visualization qualitatively matches the structure plot output from ADMIXTURE using $K=6$ (Appendix 1—figure 10). We add labels on the figure to highlight a number of pertinent features: (A) St. Lawrence Island, (B) the coastal islands and mountain ranges in British Columbia, (C) the boundary of Boreal Forest and Tundra eco-regions in the Shield Taiga, (D) Queen Elizabeth Islands, (E) Hudson Bay, and (F) Baffin Island. Many of these features were described in Schweizer et al., 2016 by interpretation of ADMIXTURE, PCA, and $F_{ST}$ statistics. FEEMS is able to succinctly provide an interpretable view of these data in a single visualization. Indeed many of these geographic features plausibly impact gray wolf dispersal and population history (Schweizer et al., 2016).

Figure 4

Download asset Open asset

We also ran EEMS on the same gray wolf dataset. We used default parameters provided by EEMS but set the number of burn-in iterations to $20\times {10}^6$, MCMC iterations to $50\times {10}^6$, and thinning intervals to 2000. We were unable to run EEMS in a reasonable run time ($\le 3$ days) for the dense spatial grid of 1207 nodes so we ran EEMS and FEEMS on a sparser graph with 307 nodes.

We find that FEEMS is multiple orders of magnitude faster than EEMS, even when running multiple runs of FEEMS for different regularization settings on both the sparse and dense graphs (Table 1). We note that constructing the graph and fitting the model with very low regularization parameters are the most computationally demanding steps in running FEEMS.

We find that many of the same geographic features that have reduced or enhanced gene-flow are concordant between the two methods. The EEMS visualization, qualitatively, best matches solutions of FEEMS with lower λ values (Figure 4, Appendix 1—figure 11); however, based on the ADMIXTURE results, visual inspection in relation to known geographical features and inspection of the observed vs fitted dissimilarity values (Appendix 1—figures 14, 22), we find these solutions to be less satisfying compared to the FEEMS solution found with λ chosen by leave-one-out cross-validation. We note that in many of the EEMS runs the MCMC appears to not have converged (based on visual inspection of trace plots) even after a large number of iterations.

See Appendix 1 ‘Mathematical notation’ for a detailed description of the notation used to describe the model. To visualize and model spatial patterns in a given population genetic dataset, FEEMS uses an undirected graph, $𝒢=(𝒱,ℰ)$ with ${\displaystyle \left|\mathcal{𝒱}\right|=d}$, where nodes represent sub-populations and edge weights ${\left(w_{k\mathrm{\ell }}\right)}_{(k,\mathrm{\ell })\in ℰ}$ represent the level of gene-flow between sub-populations $k$ and $\mathrm{\ell }$. For computational convenience, we assume $𝒢$ is a highly sparse graph, specifically a triangular grid that is embedded in geographic space around the sample coordinates. We observe a genotype matrix, $\mathit{𝒀}\in {ℝ}^{n\times p}$, with $n$ rows representing individuals and $p$ columns representing SNPs. We imagine diploid individuals are sampled on the nodes of $𝒢$ so that $y_{ij}\left(k\right)\in \{0,1,2\}$ records the count of some arbitrarily predefined allele in individual $i$, SNP $j$, on node $k\in 𝒱$. We assume a commonly used simple Binomial sampling model for the genotypes:

where conditional on $f_j\left(k\right)$ for all $j,k$, the $y_{ij}\left(k\right)$’s are independent. We then estimate an allele frequency at each node and SNP by maximum likelihood:

where $n_k$ is the number of individuals sampled at node $k$. We estimate allele frequencies at $o$ of the observed nodes out of $d$ total nodes on the graph. From Equation (1), the estimated frequency in a particular sub-population, conditional on the latent allele frequency, will approximately follow a Gaussian distribution:

Using vector notation, we represent the joint model of estimated allele frequencies as:

where ${\widehat{\mathit{𝒇}}}_j$ is a $o\times 1$ vector of estimated allele frequencies at observed nodes, ${\mathit{𝒇}}_j$ is a $d\times 1$ vector of latent allele frequencies at all the nodes (both observed and unobserved), and $\mathit{𝑨}$ is a $o\times d$ node assignment matrix where ${\mathit{𝑨}}_{k\mathrm{\ell }}=1$ if the kth estimated allele frequency comes from sub-population $\mathrm{\ell }$ and ${\mathit{𝑨}}_{k\mathrm{\ell }}=0$ otherwise; and $\text{diag}\left({\mathit{𝒅}}_{\mathit{𝒇}\mathbf{,}\mathit{𝒏}}\right)$ denotes a $o\times o$ diagonal matrix whose diagonal elements corresponds to the appropriate variance term at observed nodes.

To summarize, we estimate allele frequencies from a subset of nodes on the graph and define latent allele frequencies for all the nodes of the graph. The assignment matrix $\mathit{𝑨}$ maps these latent allele frequencies to our observations. Our summary statistics (the data) are thus $(\widehat{\mathit{𝑭}},\mathit{𝒏})$ where $\widehat{\mathit{𝑭}}$ is a $o\times p$ matrix of estimated allele frequencies and $\mathit{𝒏}$ is a $o\times 1$ vector of sample sizes for every observed node. We assume the latent allele frequencies come from a Gaussian Markov Random Field:

where $\mathit{𝑳}$ is the graph Laplacian, † represents the pseudo-inverse operator, and ${\mu }_j$ represents the average allele frequency across all of the sub-populations. Note that the multiplication by the SNP-specific factor ${\mu }_j\left(1-{\mu }_j\right)$ ensures that the variance of the latent allele frequencies vanishes as the average allele frequency approaches to 0 or 1. One interpretation of this model is that the expected squared Euclidean distance between latent allele frequencies on the graph, after being re-scaled by ${\mu }_j\left(1-{\mu }_j\right)$, is exactly the resistance distance of an electrical circuit (McRae, 2006; Hanks and Hooten, 2013):

where ${\mathit{𝒐}}_i$ is a one-hot vector (i.e. storing a 1 in element $i$ and zeros elsewhere). It is known that the resistance distance $r_{k\mathrm{\ell }}$ is equivalent to the expected commute time between nodes $k$ and $\mathrm{\ell }$ of a random walker on the weighted graph $𝒢$ (Chandra et al., 1996). Additionally, the model (Equation 3) forms a Markov random field, and thus any latent allele frequency $f_j\left(k\right)$ is conditionally independent of all other allele frequencies given its neighbors which are encoded by nonzero elements of $\mathit{𝑳}$ (Lauritzen, 1996; Koller and Friedman, 2009). Since we use a triangular grid embedded in geographic space to define the graph $𝒢$, the pattern of nonzero elements is prefixed by the structure of the sparse traingular grid.

Using the law of total variance formula, we can derive from (Equations 2, 3) an analytic form for the marginal likelihood. Before proceeding, however, we further approximate the model by assuming $\frac{1}{2}{f}_j\left(k\right)\left(1-f_j\left(k\right)\right)\approx {\sigma }^2{\mu }_j\left(1-{\mu }_j\right)$ for all $j$ and $k$ (see Appendix 1 ‘Estimating the edge weights under the exact likelihood model’ for the data model without this approximation). This assumption is mainly for computational purposes and may be a coarse approximation in general. On the other hand, the assumption is not too strong if we exclude SNPs with extremely rare allele frequencies, and more importantly, we find it leads to a good empirical performance, both statistically and computationally. With this approximation, the residual variance parameter ${\sigma }^2$ is still unknown and needs to be estimated.

Under (Equation 2, 3), the law of total variance formula leads to specific formulas for the mean and variance structure as given in (Equation 4). With those results, we arrive at the following approximate marginal likelihood:

where $\text{diag}\left({\mathit{𝒏}}^{-1}\right)$ is a $o\times o$ diagonal matrix computed from the sample sizes at observed nodes. We note the marginal distribution of ${\widehat{f}}_j$ is not necessarily a Gaussian distribution; however, we use a Gaussian approximation to facilitate computation.

To remove the SNP means we transform the estimated frequencies by a contrast matrix, $\mathit{𝑪}\in {ℝ}^{\left(o-1\right)\times o}$, that is orthogonal to the one-vector:

Let $\widehat{\mathbf{𝚺}}=\frac{1}{p}{\widehat{\mathit{𝑭}}}_{\text{s}}{\widehat{\mathit{𝑭}}}_{\text{s}}^{\top }$ be the $o\times o$ sample covariance matrix of estimated allele frequencies after re-scaling, that is, ${\widehat{\mathit{𝑭}}}_{\text{s}}$ is a matrix formed by rescaling the columns of $\widehat{\mathit{𝑭}}$ by $\sqrt{{\widehat{\mu }}_j\left(1-{\widehat{\mu }}_j\right)}$, where ${\widehat{\mu }}_j$ is an estimate of the average allele frequency (see above). We can then express the model in terms of the transformed sample covariance matrix:

where ${𝒲}_p$ denotes a Wishart distribution with $p$ degrees of freedom. Note we can equivalently use the sample squared Euclidean distance (often refereed to as a genetic distance) as a summary statistic: letting $\widehat{\mathit{𝑫}}$ be the genetic distance matrix with ${\widehat{\mathit{𝑫}}}_{k\mathrm{\ell }}={\sum }_{j=1}^p{\left({\widehat{f}}_j\left(k\right)-{\widehat{f}}_j\left(\mathrm{\ell }\right)\right)}^2/p\cdot {\widehat{\mu }}_j\left(1-{\widehat{\mu }}_j\right)$, we have

and so

using the fact that the contrast matrix $\mathit{𝑪}$ is orthogonal to the one-vector. Thus, we can use the same spatial covariance model implied by the allele frequencies once we project the distances on to the space of contrasts:

Overall, the negative log-likelihood function implied by our spatial model is the following (ignoring constant terms):

where $\mathit{𝒘}\in {ℝ}^m$ is a vectorized form of the non-zero lower-triangular entries of the weighted adjacency matrix $\mathit{𝑾}$ (recall that the graph Laplacian is completely defined by the edge weights, $\mathit{𝑳}=\text{diag}\left(\mathit{𝑾}\mathrm{𝟏}\right)-\mathit{𝑾}$, so there is an implicit dependency here). Since the graph is a triangular lattice, we only need to consider the non-zero entries to save computational time, that is, not all sub-populations are connected to each other.

We note our model (Equation 6) assumes that the $p$ SNPs are independent. This assumption is unlikely to hold when datasets are analyzed with SNPs that statistically covary (linkage disequilibrium). However, we note that the degree of freedom parameter does not affect the point estimate produced by FEEMS because it is treated as a constant term in the log-likelihood function.

One key difference between EEMS (Petkova et al., 2016) and FEEMS is how the edge weights are parameterized. In EEMS, each node is given an effective migration parameter $m_k$ for node $k\in 𝒱$ and the edge weight is parameterized as the average between the nodes it connects, that is, $w_{k\mathrm{\ell }}=\left(m_k+m_{\mathrm{\ell }}\right)/2$ for $(k,\mathrm{\ell })\in ℰ$. FEEMS, on the other hand, assigns a parameter to every nonzero edge-weight. The former has fewer parameters, with the specific consequence that it only allows isotropy and imposes an additional degree of similarity among edge weights; instead, in the latter, the edge weights are free to vary apart from the regularization imposed by the penalty. See Appendix 1 ‘Edge versus node parameterization’ and Appendix 1—figures 15, 17 for more details.

As mentioned previously, we would like to encourage that nearby edge weights on the graph have similar values to each other. This can be performed by penalizing differences between all edges connected to the same node, that is, spatially adjacent edges:

where, as before, $ℰ\left(i\right)$ denotes the set of edges that is connected to node $i$. (As mentioned earlier, in practice we choose $\alpha =1/{\widehat{w}}_0$, where ${\widehat{w}}_0$ is the solution for the ‘constant-$w$’ model, but we use the free parameter α here for full generality.) The function $x\mapsto \log \left(e^x-1\right)$ (on positive values $x\in (0,\mathrm{\infty })$) is approximately equal to $x$, for $x$ much larger than 1, and is approximately equal to $\log \left(x\right)$, for $x$ much smaller than 1. This means that our penalty function effectively penalizes differences on the log scale for edges $(i,k)$ and $(i,\mathrm{\ell })$ with very small weights, but penalizes differences on the original non-log scale for edges with large weights. Using a logarithmic-scale penalty for edges with low weights (rather than simply penalizing ${\left(w_{ik}-w_{i\mathrm{\ell }}\right)}^2$) leads to smooth graphs for small edge values, and thus allow for an additional degree of flexibility across orders of magnitude of edge weights. The penalty parameter, λ, controls the overall contribution of the penalty to the objective function. It is convenient to write the penalty in matrix-vector form which we will use throughout:

where $\mathbf{𝚫}$ is a signed graph incidence matrix derived from a unweighted graph denoting if pairs of edges are connected to the same node. Specifically, in this expression, we treat $\mathit{𝒘}$ as a vector of length $\left|ℰ\right|$ (i.e. the number of edges), and apply the function $w\mapsto \log \left(e^{\alpha w}-1\right)$ entrywise to this vector. For each pair adjacent edges $(i,k)$ and $(i,\mathrm{\ell })$ in the graph, there is a corresponding row of $\mathbf{𝚫}$ with the value +1 in the entry corresponding to edge $(i,k)$, a −1 in the entry corresponding to edge $(i,\mathrm{\ell })$, and 0’s elsewhere.

One might wonder whether it is possible to use the ${\mathrm{\ell }}_1$ norm in the penalty form Equation (8) in place of the ${\mathrm{\ell }}_2$ norm. While it is known that the ${\mathrm{\ell }}_1$ norm might increase local adaptivity and better capture the sharp changes of the underlying structure of the latent allele frequencies (e.g. Wang et al., 2016), in our case, we found an inferior performance when using the ${\mathrm{\ell }}_1$ norm over the ${\mathrm{\ell }}_2$ norm—in particular, our primary application of interest is the regime of highly missing nodes, that is, $o\ll d$, in which case the global smoothing seems somewhat necessary to encourage stable recovery of the edge weights at regions with sparsely observed nodes (see Appendix 1 ‘Smooth penalty with ${\mathrm{\ell }}_{\mathrm{1}}$norm’). In addition, adding the penalty ${\varphi }_{\lambda ,\alpha }\left(\mathit{𝒘}\right)$ allows us to implement faster algorithms to solve the optimization problem due to the differentiability of the ${\mathrm{\ell }}_2$ norm, and as a result, it leads to better overall computational savings and a simpler implementation.

Putting Equation (7) and Equation (8) together, we infer the migration edge weights $\widehat{\mathit{𝒘}}$ by minimizing the following penalized negative log-likelihood function:

where $\mathit{𝒍},\mathit{𝒖}\in {ℝ}_{+}^m$ represent respectively the entrywise lower- and upper bounds on $\mathit{𝒘}$, that is, we constrain the lower- and upper bound of the edge weights to $\mathit{𝒍}$ and $\mathit{𝒖}$ throughout the optimization. When no prior information is available on the range of the edge weights, we often set $\mathit{𝒍}=\mathrm{}$ and $\mathit{𝒖}=+\mathrm{\infty }$.

One advantage of the formulation of Equation 9 is the use of the vector form parameterization $\mathit{𝒘}\in {ℝ}_{+}^m$ of the symmetric weighted adjacency matrix $\mathit{𝑾}\in {ℝ}_{+}^{d\times d}$. In our triangular graph $𝒢=(𝒱,ℰ)$, the number of non-zero lower-triangular entries is $m=𝒪\left(d\right)\ll d^2$, so working directly on the space of vector parameterization saves computational cost. In addition, this avoids the symmetry constraint imposed on the adjacency matrix $\mathit{𝑾}$, hence making optimization easier (Kalofolias, 2016).

We solve the optimization problem using a constrained quasi-Newton optimization algorithm, specifically L-BFGS implemented in scipy (Byrd et al., 1995; Virtanen et al., 2020) (RRID:SCR_008058). Since our objective Equation 9 is non-convex, the L-BFGS algorithm is guaranteed to converge only to a local minimum. Even so, we empirically observe that local minima starting from different initial points are qualitatively similar to each other across many datasets. The L-BFGS algorithm requires gradient and objective values as inputs. Note the naive computation of the objective Equation 9 is computationally prohibitive because inverting the graph Laplacian has complexity $𝒪\left(d^3\right)$. We take advantage of the sparsity of the graph and specific structure of the problem to efficiently compute gradient and objective values. In theory, our implementation has computational complexity of $𝒪\left(do+o^3\right)$ per iteration which, in the setting of $o\ll d$, is substantially smaller than $𝒪\left(d^3\right)$. It is possible to achieve $𝒪\left(do+o^3\right)$ per-iteration complexity by using a solver that is specially designed for a sparse Laplacian system. In our work, we use sparse Cholesky factorization which may slightly slow down the per-iteration complexity (See Appendix Material for the details of the gradient and objective computation).

For estimating the residual variance parameter ${\sigma }^2$, we first estimate it via maximum likelihood assuming homogeneous isolation by distance. This corresponds to the scenario where every edge-weight in the graph is given the exact same unknown parameter value $w_0$. Under this model, we only have two unknown parameters $w_0$ and the residual variance ${\sigma }^2$. We estimate these two parameters by jointly optimizing the marginal likelihood using a Nelder-Mead algorithm implemented in scipy (Virtanen et al., 2020) (RRID:SCR_008058). This requires only likelihood computations which are efficient due to the sparse nature of the graph. This optimization routine outputs an estimate of the residual variance ${\widehat{\sigma }}^2$ and the null edge weight ${\widehat{w}}_0$, which can be used to construct $\mathit{𝑾}\left({\widehat{w}}_0\right)$ and in turn $\mathit{𝑳}\left({\widehat{w}}_0\right)$.

One strategy we found effective is to fit the model of homogeneous isolation by distance and then fix the estimated residual variance ${\widehat{\sigma }}^2$ throughout later fits of the more flexible penalized models—See Appendix 1 ‘Jointly estimating the residual variance and edge weights’. Additionally, we find that initializing the edge weights to ${\widehat{w}}_0$ to be a useful and intuitive strategy to set the initial values for the entries of $\mathit{𝒘}$ to the correct scale.

FEEMS estimates one set of graph edge weights for each setting of the tuning parameters λ and α which control the smoothness of the fitted edge weights. Figure 3 shows that the estimated migration surfaces vary substantially depending on the particular choices of the tuning parameters, and indeed, due to the large fraction of unobserved nodes, it can highly over-fit the observed data unless regularized accordingly. To address the issue of selecting the tuning parameters, we propose using leave-one-out cross-validation to assess each fitted model’s generalization ability at held out locations.

To simplify the notation, we write the model Equation 4 for the estimated allele frequencies in SNP $j$ as

where

For each fold, we hold out one node from the set of observed nodes in the graph and use the rest of the nodes to fit FEEMS across a sequential grid of regularization parameters. Note that our objective function is non-convex, so the algorithm converges to different local minima for different regularization parameters, even with the same initial value ${\widehat{w}}_0$. To stabilize the cross-validation procedure, we recommend using a warm start strategy in which one solves the problem for the largest value of regularization parameters first and use this solution to initialize the algorithm at the next largest value of regularization parameters, and so on. Empirically, we find that using warm starts gives far more reliable model selection than with cold starts, where the problems over the sequence of parameters are solved independently with same initial value ${\widehat{w}}_0$. We suspect that the poor performance of leave-one-out cross-validation without warm starts is attributed to spatial dependency of allele frequencies and the large fraction of unobserved nodes. Without loss of generality, we assume that the last node has been held out. Re-writing the distribution of the observed frequencies according to the split of observed nodes,

the conditional mean of the observed frequency ${\widehat{f}}_j^{\text{val}}$ on the held out node, given the rest, is given by

Using this formula, we can predict allele frequencies at held out locations using the fitted graph $\widehat{\mathit{𝑳}}=\widehat{\mathit{𝑳}}(\lambda ,\alpha )$ for each setting of tuning parameters λ and α. Note that in Equation (10), the parameters ${\mu }_j$ and σ are also unknown, and we use an estimate of the average allele frequency ${\widehat{\mu }}_j$ and the estimated residual variance $\widehat{\sigma }$ from the ‘constant-$w$’ model (they are not dependent on λ and α). Then we select the tuning parameters λ and α that output the minimum prediction error averaged over all SNPs ${\displaystyle \frac{1}{p}\sum _j{\Vert {\hat{f}}_j^{\text{val,pred}}-{\hat{f}}_j^{\text{val}}\Vert }_2^2}$, averaged over all the held out nodes (with $o$ observed nodes in total). As mentioned earlier, in practice we choose $\alpha =1/{\widehat{w}}_0$ and hence we can use the leave-one-out cross-validation to search for λ only, which allows us to avoid the computational cost of searching over the two-dimensional parameter space.

At a high level, we can summarize the differences between FEEMS and EEMS as follows: (1) the likelihood functions of FEEMS and EEMS are slightly different as a function of the graph Laplacian $\mathit{𝑳}$; (2) the migration rates are parameterized in terms of edge weights or in terms of node weights; and (3) EEMS is based on Bayesian inference and thus chooses a prior and studies the posterior distribution, while FEEMS is an optimization-based approach and thus chooses a penalty function and minimizes the penalized log-likelihood (in particular, the EEMS prior and the FEEMS penalty are both aiming for locally constant type migration surfaces). The last two points were already discussed in the above sections, so here we focus on the difference of the likelihoods between the two methods.

FEEMS develops the spatial model for the genetic differentiation through Gaussian Markov Random Field, but the resulting likelihood has similarities to EEMS (Petkova et al., 2016) which considers the pairwise coalescent times. Using our notation, we can write the EEMS model as

where $\nu \in [o-1,p]$ is the effective degree of freedom, ${\sigma }^{*}>0$ is the scale nuisance parameter, and $\mathit{𝒒}$ is a $d\times 1$ vector of the within-sub-population coalescent rates. ${\widehat{\mathit{𝑫}}}^{*}$ represents the genetic distance matrix without re-scaling, where the $(k,\mathrm{\ell })$-th element is given by ${\widehat{\mathit{𝑫}}}_{k\mathrm{\ell }}^{\star }={\sum }_{j=1}^p{\left({\widehat{f}}_j\left(k\right)-{\widehat{f}}_j\left(\mathrm{\ell }\right)\right)}^2/p$. That is, unlike FEEMS, EEMS does not consider the SNP-specific re-scaling factor ${\mu }_j\left(1-{\mu }_j\right)$ to account for the vanishing variance of the observed allele frequencies as the average allele frequency approaches to 0 or 1.

In Equation (11), the effective degree of freedom ν is introduced to account for the dependency across SNPs in close proximity. Because EEMS uses a hierarchical Bayesian model to infer the effective migration rates, ν is being estimated alongside other model parameters. On the other hand, FEEMS uses an optimization-based approach and the degrees of freedom has no influence on the point estimate of the migration rates. Besides the effective degree of freedom and the SNP-specific re-scaling by ${\mu }_j\left(1-{\mu }_j\right)$, the EEMS and FEEMS likelihoods are equivalent up to constant factors, as long as only one individual is observed per node and the residual variance ${\sigma }^2$ is allowed to vary across nodes—See Appendix 1 ‘Jointly estimating the residual variance and edge weights’ for details. The constant factors, such as ${\sigma }^{*}$, can be effectively absorbed into the unknown model parameters $\mathit{𝑳}$ and $\mathit{𝒒}$ and therefore they do not affect the estimation of effective migration rates, up to constant factors.

We analyzed a population genetic dataset of North American gray wolves previously published in Schweizer et al., 2016. For this, we downloaded plink (RRID:SCR_001757) formatted files and spatial coordinates from https://doi.org/10.5061/dryad.c9b25. We removed all SNPs with minor allele frequency less than 5% and with missingness greater then 10%, resulting in a final set of 111 individuals and 17,729 SNPs.

We fit the Pritchard, Donnelly, and Stephens model (PSD) and ran principal components analysis on the genotype matrix of North American gray wolves (Price et al., 2006; Pritchard et al., 2000). For the PSD model, we used the ADMIXTURE software (RRID:SCR_001263) on the un-normalized genotypes, running five replicates per choice of $K$, from $K=2$ to $K=8$ (Alexander et al., 2009). For each $K$, we choose the one that achieved the highest likelihood to visualize. For PCA, we centered and scaled the genotype matrix and then ran sklearn (RRID:SCR_019053) implementation of PCA, truncated to compute 50 eigenvectors.

To create a dense triangular lattice around the sample locations, we first define an outer boundary polygon. As a default, we construct the lattice by creating a convex hull around the sample points and manually trimming the polygon to adhere to the geography of the study organism and balancing the sample point range with the extent of local geography using the following website https://www.keene.edu/campus/maps/tool/. We often do not exclude internal ‘holes’ in the habitat (e.g. water features for terrestrial animals), and let the model instead fit effective migration rates for those features to the extent they lead to elevated differentiation. We also emphasize the importance of defining the lattice for FEEMS as well as EEMS and suggest this should be carefully curated with prior biological knowledge about the system.

To ensure edges cover an equal area over the entire region, we downloaded and intersected a uniform grid defined on the spherical shape of earth (Sahr et al., 2003). These defined grids are pre-computed at a number of different resolutions, allowing a user to test FEEMS at different grid densities which is an important feature to explore.

The code to reproduce the results of this paper and more can be found at https://github.com/jhmarcus/feems-analysis (Marcus and Ha, 2021a, copy archived at swh:1:rev:f2d7330f25f8a11124db09000918ae38ae00d4a7, Marcus and Ha, 2021b). A python (RRID:SCR_008394) package implementing the method can be found at https://github.com/Novembrelab/feems.

We denote matrices using bold capital letters $\mathit{𝑨}$. Bold lowercase letters are vectors $\mathit{𝒂}$, and non-bold lowercase letters are scalars $a$. We denote by ${\mathit{𝑨}}^{-1}$ and ${\mathit{𝑨}}^{†}$ the inverse and (Moore-Penrose) pseudo-inverse of $\mathit{𝑨}$, respectively. We use $\mathit{𝒚}\sim N_p(\mathit{𝝁},\mathbf{𝚺})$ to express that the random vector $\mathit{𝒚}$ is modeled as a $p$-dimensional multivariate Gaussian distribution with fixed parameters $\mathit{𝝁}$ and $\mathbf{𝚺}$ and use the conditional notation $\mathit{𝒚}|\mathit{𝝁}\sim N_p(\mathit{𝝁},\mathbf{𝚺})$ if $\mathit{𝝁}$ is random.

A graph is a pair $𝒢=(𝒱,ℰ)$, where $𝒱$ denotes a set of nodes or vertices and $ℰ\subseteq 𝒱\times 𝒱$ denotes a set of edges. Throughout we assume the graph $𝒢$ is undirected, weighted, and contains no self loops, that is, $(k,\mathrm{\ell })\in ℰ\iff (\mathrm{\ell },k)\in ℰ$ and $(k,k)\notin ℰ$ and each edge $(k,\mathrm{\ell })\in ℰ$ is given a weight $w_{k\mathrm{\ell }}=w_{\mathrm{\ell }k}>0$. We write $\mathit{𝑾}$ to indicate the symmetric weighted adjacency matrix, that is,

$\mathit{𝒘}\in {ℝ}^m$ is a vectorized form of the non-zero lower-triangular entries of $\mathit{𝑾}$ where ${\displaystyle m=\left|\mathcal{ℰ}\right|/2}$ is the number of non-zero lower triangular elements. We denote by $\mathit{𝑳}=\text{diag}\left(\mathit{𝑾}\mathrm{𝟏}\right)-\mathit{𝑾}$ the graph Laplacian.

In practice, we make a change of variable from $\mathit{𝒘}\in {ℝ}_{+}^m$ to $\mathit{𝒛}=\log \left(\mathit{𝒘}\right)\in {ℝ}^m$ and the algorithm is applied to the transformed objective function:

After the change of variable, the objective value remains the same, whereas it follows from the chain rule that $\nabla \left(\tilde{\mathrm{\ell }}\left(\mathit{𝒛}\right)+{\tilde{\varphi }}_{\lambda ,\alpha }\left(\mathit{𝒛}\right)\right)=\nabla \left(\mathrm{\ell }\left(\mathit{𝒘}\right)+{\varphi }_{\lambda ,\alpha }\left(\mathit{𝒘}\right)\right)\odot \mathit{𝒘}$ where $\odot$ indicates the Hadamard product or elementwise product—for notational convenience, we drop the dependency of $\mathrm{\ell }$ on the quantities ${\sigma }^2$ and $\mathit{𝑪}\widehat{\mathbf{𝚺}}{\mathit{𝑪}}^{\top }$. Furthermore, the computation of $\nabla {\varphi }_{\lambda ,\alpha }\left(\mathit{𝒘}\right)$ is relatively straightforward, so in the rest of this section, we discuss only the computation of the gradient of the negative log-likelihood function with respect to $\mathit{𝒘}$, that is, $\nabla \mathrm{\ell }\left(\mathit{𝒘}\right)$.

Recall, by definition, the graph Laplacian $\mathit{𝑳}$ implicitly depends on the variable $\mathit{𝒘}$ through $\mathit{𝑳}=\text{diag}\left(\mathit{𝑾}\mathrm{𝟏}\right)-\mathit{𝑾}$. Throughout we assume the first $o$ rows and columns of $\mathit{𝑳}$ correspond to the observed nodes. With this assumption, our node assignment matrix has block structure $\mathit{𝑨}=\left[{\mathbf{𝐈}}_{o\times o}|{\mathbf{\hspace{0.33em}0}}_{o\times \left(d-o\right)}\right]$. To simplify some of the equations appearing later, we introduce the notation: we define

and

Applying the chain rule and matrix derivatives, we can calculate:

where vec is the vectorization operator and $\partial \mathrm{\ell }/\partial \text{vec}\left(\mathit{𝑳}\right)$ and $\partial \text{vec}\left(\mathit{𝑳}\right)/\partial {\mathit{𝒘}}^{\top }$ are $1\times d^2$ vector and $d^2\times d$ matrix, respectively, given by

Here, $\mathit{𝑺}$ and $\mathit{𝑻}$ are linear operators that satisfy $\mathit{𝑺}\mathit{𝒘}=\text{diag}\left(\mathit{𝑾}\mathrm{𝟏}\right)$ and $\mathit{𝑻}\mathit{𝒘}=\mathit{𝑾}$. Note $\mathit{𝑺}$ and $\mathit{𝑻}$ both have $𝒪\left(d\right)$ many nonzero entries, so we can perform sparse matrix multiplication to efficiently compute the matrix-vector multiplication $\partial \mathrm{\ell }/\partial \text{vec}\left(\mathit{𝑳}\right)\cdot \left(\mathit{𝑺}-\mathit{𝑻}\right)$. On the other hand, the computation of $\partial \mathrm{\ell }/\partial \text{vec}\left(\mathit{𝑳}\right)$ is more challenging as it requires inverting the full $d\times d$ matrix ${\mathit{𝑳}}_{\text{full}}$. Next, we develop a procedure that efficiently computes $\partial \mathrm{\ell }/\partial \text{vec}\left(\mathit{𝑳}\right)$. We proceed by dividing the task into multiple steps.

The graph Laplacian $\mathit{𝑳}$ is orthogonal to the one vector 1, so using the notation introduced in Equation (12), we can express our objective function as

With the identity $\mathbf{𝐈}-{\mathbf{𝚷}}_{\mathrm{𝟏}}=\mathbf{𝚺}{\mathit{𝑪}}^{\top }{\left(\mathit{𝑪}\mathbf{𝚺}{\mathit{𝑪}}^{\top }\right)}^{-1}\mathit{𝑪}$, the trace term is:

The matrix inside the trace has been constructed in the gradient computation, see Equation (18). In terms of the determinant, we use the same approach considered in Petkova et al., 2016—in particular, concatenating ${\mathit{𝑪}}^{\top }$ and $\mathrm{𝟏}$, the matrix $\left[{\mathit{𝑪}}^{\top }|\mathbf{\hspace{0.33em}1}\right]$ is orthogonal, so it can be shown that

Rearranging terms and using the fact $det\left({\mathit{𝑼}}^{-1}\right)=det{\left(\mathit{𝑼}\right)}^{-1}$ for any matrix $\mathit{𝑼}$, we obtain:

We have computed ${\mathbf{𝚺}}^{-1}$ in Equation (14), so each of the terms above can be computed without any additional matrix multiplications. Finally, the signed graph incidence matrix $\mathbf{𝚫}$ defined on the edges of the graph is, by construction, highly sparse with $𝒪\left(d\right)$ many nonzero entries. Hence we implement sparse matrix multiplication to evaluate the penalty function ${\varphi }_{\lambda ,\alpha }\left(\mathit{𝒘}\right)$ while avoiding the full-dimensional matrix-vector product.

When we developed the FEEMS model, we used the approximation $\frac{1}{2}{f}_j\left(k\right)\left(1-f_j\left(k\right)\right)\approx {\sigma }^2{\mu }_j\left(1-{\mu }_j\right)$ for all SNPs $j$ and all nodes $k$ (see Equation 4) and estimated the residual variance ${\sigma }^2$ under the homogeneous isolation by distance model. The primary reason of using this approximation was primarily computational. While the approximation is not too strong if SNPs with rare allele frequencies are excluded, it is also critical that the estimation quality of the migration rates is not affected. In this subsection we introduce the inferring procedure of the migration rates under the exact likellihood model and compare it with FEEMS.

Note that without approximation, we can calculate the exact analytical form for the marginal likelihood of the estimated frequency as follows (after removing the SNP means):

where ${\left\{a_k\right\}}_{k=1}^d$ represents the vector $\mathit{𝒂}=(a_1,\mathrm{\dots },a_d)$. Compared to the model Equation (5), this expression does not introduce the unknown residual variance parameter ${\sigma }^2$ and instead each node has its own residual parameter given by $\left(1-L_{kk}^{†}\right)/2$. Because the residual parameters must be positive, this means that we have to search for the graphs that ensure $L_{kk}^{†}\le 1$ for all nodes $k$. With that said, we can consider the following constrained optimization problem:

where ${\mathrm{\ell }}_{\text{𝖾𝗑𝖺𝖼𝗍}}$ is the negative log-likelihood function based on Equation (19) and ${\varphi }_{\lambda ,\alpha }$ is the smooth penalty function defined earlier. The main difficulty of solving Equation (20) is that enforcing the constraint $L_{kk}^{†}\le 1$ for all nodes $k\in 𝒱$, requires full computation of the pseudo-inverse of a $d\times d$ matrix $\mathit{𝑳}$ which is computationally demanding. We instead relax the constraint and consider the following form as a proxy for optimization Equation (20):

Note that the constraint $L_{kk}^{†}\le 1$ is now placed at the observed nodes only, which can lead to computational savings if $o\ll d$. The problem Equation (21) can be solved efficiently using any gradient-based algorithms where we can calculate the gradient of ${\mathrm{\ell }}_{\text{𝖾𝗑𝖺𝖼𝗍}}$ with respect to $\mathit{𝑳}$ as

where $\mathit{𝑴}$ is a $o\times o$ matrix defined in Equation (18), and $\mathit{𝑵}$ is a $o\times d^2$ matrix whose rows correspond to the observed subsets of the rows of the $d^2\times d^2$ matrix ${\mathit{𝑳}}_{\text{full}}^{-1}\otimes {\mathit{𝑳}}_{\text{full}}^{-1}$.

Appendix 1—figure 12 shows the result when the penalized maximum likelihood Equation (21) is applied to the North American wolf dataset with a setting of $\lambda =2.06$ (the same value of λ as given in Figure 4) and $\alpha =1/{\widehat{w}}_0$, where ${\widehat{w}}_0$ is the solution for the ‘constant-$w$’ model. We can see that the resulting estimated migration surfaces are qualitatively similar to that shown in Figure 4. We also observed similar results between FEEMS and the penalized maximum likelihood Equation (21) across multiple datasets. On the other hand, we found that at the fitted surface the residual variances ${\displaystyle 1-L_{kk}^{†}}$ are not always positive because the constraints are enforced only at the observed nodes. This is problematic because it can cause the model to be ill-defined at the unobserved nodes and make the algorithm numerically unstable. Note that FEEMS avoids this issue by decoupling the residual variance parameter ${\sigma }^2$ from the graph-related parameters ${\displaystyle \mathit{w}}$. The resulting model Equation (6) also has more resemblance to spatial coalescent model used in EEMS (Petkova et al., 2016), and we thus recommend using FEEMS as a primary method for inferring migration rates.

One simple strategy we have used throughout the paper was to fit ${\sigma }^2$ first under a model of homogeneous isolation by distance and prefix the estimated residual variance to the resulting ${\widehat{\sigma }}^2$ for later fits of the effective migration rates. Alternatively, we can consider estimating the unknown residual variance simultaneously with the edge weights, instead of prefixing it from the estimation of the null model—the hope here is to simultaneously correct the model misspecification and allow for improving model fit to the data. To develop the framework for simultaneous estimation of the residual variance and edge weights, let us consider a model that generalizes both Equation (6) and Equation (19), that is,

where ${\mathit{𝝈}}^2$ is a $d\times 1$ vector of node specific residual variance parameters, that is, each deme has its own residual parameter ${\sigma }_k$. If the parameters ${\sigma }_k$’s are assumed to be the same across nodes, this reduces to the FEEMS model Equation (6) while setting ${\displaystyle {\sigma }_k=(1-{\mathit{L}}_{kk}^{†})/2}$ gives the model Equation (19). Then we solve the following optimization problem

where ${\mathrm{\ell }}_{\text{𝗃𝗈𝗂𝗇𝗍}}$ is the negative log-likelihood function based on Equation (22). Note that the residual variances and edge weights are both searched in the optimization for finding the optimal solutions. To solve the problem, we can use the quasi-newton algorithm for optimizing the objective function.

Appendix 1—figure 13 shows the fitted graphs with different strategies of estimating the residual variances. Appendix 1—figure 13A shows the result when the model has a single residual variance ${\sigma }^2$, and Appendix 1—figure 13B shows the result when the residual variances are allowed to vary across nodes. In both cases, estimating the residual variances jointly with the edge weights yields similar and comparable outputs to the default setting of prefixing it from the null model (Figure 4), except that we can further observe reduced effective migration around Queen Elizabeth Islands as shown in Appendix 1—figure 13B. In EEMS, in order to estimate the genetic diversity parameters for every spatial location, which play a similar role as the residual variances in FEEMS, a Voronoi-tessellation prior is placed to encourage sharing of information across adjacent nodes and prevent over-fitting. Similarly, we can place the spatial smooth penalty on the residual variances (i.e. ${\displaystyle {\varphi }_{\lambda ,\alpha }}$ defined on the variable ${\mathit{𝝈}}^2)$, but it introduces additional hyperparameters to tune, without substantially improving the model’s fit to the data. In this work, we choose to fit the single residual variance ${\sigma }^2$ under the null model and prefix it as a simple but effective strategy with apparent good empirical performance.

One of the novel features of FEEMS is its ability to directly fit the edge weights of the graph that best suit the data. This direct edge parameterization may increase the risk of model’s overfitting, but also allows for more flexible estimation of migration histories. Furthermore, as seen in Figure 2 and Appendix 1—figure 2, it has potential to recover anisotropic migration processes. This is in contrast to EEMS wherein every spatial node is assigned an effective migration parameter $m_k$ and a migration rate on each edge joining nodes $k$ and $\mathrm{\ell }$ is given by the average effective migration ${\displaystyle w_{k\ell }=(m_k+m_{\ell })/2}$. Not surprisingly, by assigning each edge to be the average of connected nodes, a form of implicit spatial regularization is imposed because multiple edges connected to the same node would average that node’s parameter value. In some cases, this has the desirable property of imposing an additional degree of similarity across edge weights, but at the same time it also restricts the model’s capacity to capture a richer set of structure present in the data (e.g. Petkova et al., 2016, Supplementary Figure 2). To be concrete, Appendix 1—figure 15 displays two different fits of FEEMS based on edge parameterization (Appendix 1—figure 15A) and node parameterization (Appendix 1—figure 15B), run on a previously published dataset of human genetic variation from Africa (see Peter et al., 2020 for details on the description of the dataset). Running FEEMS with a node-based parameterization is straightforward in our framework—all we have to do is to reparameterize the edge weights by the average effective migration and solve the corresponding optimization problem (Optimization) with respect to $\mathit{𝒎}$. It is evident from the results that FEEMS with edge parameterization exhibits subtle correlations that exist between the annotated demes in the figure, whereas node parameterization fails to recover them. We also compare the model fit of FEEMS to the observed genetic distance (Appendix 1—figure 16) and find that edge-based parameterization provides a better fit to the African dataset. Appendix 1—figure 17 further demonstrates that in the coalescent simulations with anisotropic migration, the node parameterization is unable to recover the ground truth of the underlying migration rates even when the nodes are fully observed.

FEEMS’s primary optimization objective (see Equation 9) is:

where the spatial smoothness penalty is given by an ${\mathrm{\ell }}_2$-based penalty function: ${\displaystyle {\varphi }_{\lambda ,\alpha }(\mathit{w})=\frac{\lambda }{2}{\Vert \mathbf{\Delta }\log (e^{\alpha \mathit{w}}-\mathbf{1})\Vert }_2^2}$. It is well known that an ${\mathrm{\ell }}_1$-based penalty can lead to a better local adaptive fitting and structural recovery than ${\mathrm{\ell }}_2$-based penaltyies (Wang et al., 2016), but at the cost of handling non-smooth objective functions that are often computationally more challenging. In a spatial genetic dataset, one major challenge is to deal with the relatively sparse sampling design where there are many unobserved nodes on the graph. In this statistically challenging scenario, we found that an ${\mathrm{\ell }}_2$-based penalty allows for more accurate and reliable estimation of the geographic features.

Specifically, writing ${\displaystyle {\varphi }_{\lambda ,\alpha }^{{\ell }_1}(\mathit{w})=\lambda {\Vert \mathbf{\Delta }\log (e^{\alpha \mathit{w}}-\mathbf{1})\Vert }_1}$, we considered the alternate following composite objective function: