Selective Inference for Testing Trees and Edges in Phylogenetics

Selective inference is considered for testing trees and edges in phylogenetic tree selection from molecular sequences. This improves the previously proposed approximately unbiased test by adjusting the selection bias when testing many trees and edges at the same time. The newly proposed selective inference $p$-value is useful for testing selected edges to claim that they are significantly supported if $p>1-\alpha$, whereas the non-selective $p$-value is still useful for testing candidate trees to claim that they are rejected if $p<\alpha$. The selective $p$-value controls the type-I error conditioned on the selection event, whereas the non-selective $p$-value controls it unconditionally. The selective and non-selective approximately unbiased $p$-values are computed from two geometric quantities called signed distance and mean curvature of the region representing tree or edge of interest in the space of probability distributions. These two geometric quantities are estimated by fitting a scaling-law model to the non-parametric multiscale bootstrap probabilities. For better understanding the geometric nature of the problem, a visualization of probability distributions is presented. Our general method is applicable to a wider class of problems; phylogenetic tree selection is an example of model selection, and it is interpreted as the variable selection of multiple regression, where each edge corresponds to each predictor. Our method is illustrated in a previously controversial phylogenetic analysis of human, rabbit and mouse.


INTRODUCTION
A phylogenetic tree is a diagram showing evolutionary relationships among species, and a tree topology is a graph obtained from the phylogentic tree by ignoring the branch lengths. The primary objective of any phylogenetic analysis is to approximate a topology that reflects the evolution history of the group of organisms under study. Branches of the tree are also referred to as edges in the tree topology. Given a rooted tree topology, or a unrooted tree topology with an outgroup, each edge splits the tree so that it defines the clade consisting of all the descendant species. Therefore, edges in a tree topology represent clades of species. Because the phylogenetic tree is commonly inferred from molecular sequences, it is crucial to assess the statistical confidence of the inference. In phylogenetics, it is a common practice to compute confidence levels for tree topologies and edges. For example, the bootstrap probability (Felsenstein, 1985) is the most commonly used confidence measure, and other methods such as the Shimodaira-Hasegawa test (Shimodaira and Hasegawa, 1999) and the multiscale bootstrap method (Shimodaira, 2002) are also often used. However, these conventional methods are limited in how well they address the issue of multiplicity when there are many alternative topologies and edges. Herein, we discuss a new approach, selective inference (SI), that is designed to address the issue of multiplicity.
For illustrating the idea of selective inference, we first look at a simple example of 1-dimensional normal random variable Z with unknown mean θ ∈ R and variance 1: Z ∼ N(θ , 1). (1) Observing Z = z, we would like to test the null hypothesis H 0 : θ ≤ 0 against the alternative hypothesis H 1 : θ > 0. We denote the cumulative distribution function of N(0, 1) as (x) and define the upper tail probability as¯ (x) = 1 − (x) = (−x). Then, the ordinary (i.e., non-selective) inference leads to the p-value of the one-tailed z-test as p(z) : = P(Z > z | θ = 0) =¯ (z). (2) What happens when we test many hypotheses at the same time? Consider random variables Z i ∼ N(θ i , 1), i = 1, . . . , K all , not necessarily independent, with null hypotheses θ i ≤ 0, where K true hypotheses are actually true. To control the number of falsely rejecting the K true hypotheses, there are several multiplicity adjusted approaches such as the family-wise error rate (FWER) and the false discovery rate (FDR). Instead of testing all the K all hypotheses, selective inference (SI) allows for K select hypotheses with z i > c i for constants c i specified in advance. This kind of selection is very common in practice (e.g., publication bias), and it is called as the file drawer problem by Rosenthal (1979). Instead of controlling the multiplicity of testing, SI alleviates it by reducing the number of tests. The mathematical formulation of SI is easier than FWER and FDR in the sense that hypotheses can be considered separately instead of simultaneously. Therefore, we simply write z > c by dropping the index i for one of the hypotheses. In selective inference, the selection bias is adjusted by considering the conditional probability given the selection event, which leads to the following p-value (Fithian et al., 2014;Tian and Taylor, 2018) p(z, c) : = P(Z > z | Z > c, θ = 0) =¯ (z)/¯ (c), where p(z) of Equation (2) is divided by the selection probability P(Z > c | θ = 0) =¯ (c). In the case of c = 0, this corresponds to the two-tailed z-test, because the selection probability is¯ (0) = 0.5 and p(z, c) = 2p(z). For significance level α (we use α = 0.05 unless otherwise stated), it properly controls the type-I error conditioned on the selection event as P(p(Z, c) < α | Z > c, θ = 0) = α, while the non-selective p-value violates the type-I error as P(p(Z) < α | Z > c, θ = 0) = α/¯ (c) > α. The selection bias can be very large when¯ (c) ≪ 1 (i.e., c ≫ 0), or K select ≪ K all . Selective inference has been mostly developed for inferences after model selection (Taylor and Tibshirani, 2015;Tibshirani et al., 2016), particularly variable selection in regression settings such as lasso (Tibshirani, 1996). Recently, Terada and Shimodaira (2017) developed a general method for selective inference by adjusting the selection bias in the approximately unbiased (AU) p-value computed by the multiscale bootstrap method (Shimodaira, 2002(Shimodaira, , 2004(Shimodaira, , 2008. This new method can be used to compute, for example, confidence intervals of regression coefficients in lasso (Figure 1). In this paper, we apply this method to phylogenetic inference for computing proper confidence levels of tree topologies (dendrograms) and edges (clades or clusters) of species. As far as we know, this is the first attempt to consider selective inference in phylogenetics. Our selective inference method is implemented in software scaleboot (Shimodaira, 2019) working jointly with CONSEL (Shimodaira and Hasegawa, 2001) for phylogenetics, and it is also implemented in a new version of pvclust (Suzuki and Shimodaira, 2006) for hierarchical clustering, where only edges appeared in the observed tree are "selected" for computing p-values. Although our argument is based on the rigorous theory of mathematical statistics in Terada and Shimodaira (2017), a self-contained illustration is presented in this paper for the theory as well as the algorithm of selective inference.
Phylogenetic tree selection is an example of model selection. Since each tree can be specified as a combination of edges, tree selection can be interpreted as the variable selection of multiple regression, where edges correspond to the predictors of regression (Shimodaira, 2001;Shimodaira and Hasegawa, 2005). Because all candidate trees have the same number of model parameters, the maximum likelihood (ML) tree is obtained by comparing log-likelihood values of trees (Felsenstein, 1981). In order to adjust the model complexity by the number of parameters in general model selection, we compare Akaike Information Criterion (AIC) values of candidate models (Akaike, 1974). AIC is used in phylogenetics for selecting the substitution model (Posada and Buckley, 2004). There are several modifications of AIC that allow for model selection. These include the precise estimation of the complexity term known as Takeuchi Information Criterion (Burnham and Anderson, 2002;Konishi and Kitagawa, 2008), and adaptations for incomplete data (Shimodaira and Maeda, 2018) and covariate-shift data (Shimodaira, 2000). AIC and all these modifications are derived FIGURE 2 | Examples of two unrooted trees T1 and T7. Branch lengths represent ML estimates of parameters (expected number of substitutions per site). T1 includes edges E1, E2, and E3 and T7 includes E1, E6, and E8.
for estimating the expected Kullback-Leibler divergence between the unknown true distribution and the estimated probability distribution on the premise that the model is misspecified. When using regression model for prediction purpose, it may be sufficient to find only the best model which minimizes the AIC value. Considering random variations of dataset, however, it is obvious in phylogenetics that the ML tree does not necessarily represent the true history of evolution. Therefore, Kishino and Hasegawa (1989) proposed a statistical test whether two log-likelihood values differ significantly (also known as Kishino-Hasegawa test). The log-likelihood difference is often not significant, because its variance can be very large for non-nested models when the divergence between two probability distributions is large; see Equation (26) in section 6.1. The same idea of model selection test whether two AIC values differ significantly has been proposed independently in statistics (Linhart, 1988) and econometrics (Vuong, 1989). Another method of model selection test (Efron, 1984) allows for the comparison of two regression models with an adjusted bootstrap confidence interval corresponding to the AU p-value. For testing which model is better than the other, the null hypothesis in the model selection test is that the two models are equally good in terms of the expected value of AIC on the premise that both models are misspecified. Note that the null hypothesis is whether the model is correctly specified or not in the traditional hypothesis testing methods including the likelihood ratio test for nested models and the modified likelihood ratio test for non-nested models (Cox, 1962). The model selection test is very different from these traditional settings. For comparing AIC values of more than two models, a multiple comparisons method is introduced to the model selection test (Shimodaira, 1998;Shimodaira and Hasegawa, 1999), which computes the confidence set of models. But the multiple comparisons method is conservative by nature, leading to more false negatives than expected, because it considers the worst scenario, called the least favorable configuration. On the other hand, the model selection test (designed for two models) and bootstrap probability (Felsenstein, 1985) lead to more false positives than expected when comparing more than two models (Shimodaira and Hasegawa, 1999;Shimodaira, 2002). The AU p-value mentioned earlier has been developed for solving this problem, and we are going to upgrade it for selective inference.

PHYLOGENETIC INFERENCE
For illustrating phylogenetic inference methods, we analyze a dataset consisting of mitochondrial protein sequences of six mammalian species with n = 3, 414 amino acids (n is treated as sample size). The taxa are labeled as 1=Homo sapiens (human), 2=Phoca vitulina (seal), 3=Bos taurus (cow), 4=Oryctolagus cuniculus (rabbit), 5=Mus musculus (mouse), and 6=Didelphis virginiana (opossum). The dataset will be denoted as X n = (x 1 , . . . , x n ). The software package PAML (Yang, 1997) was used to calculate the site-wise log-likelihoods for trees. The mtREV model (Adachi and Hasegawa, 1996) was used for amino acid substitutions, and the site-heterogeneity was modeled by the discrete-gamma distribution (Yang, 1996). The dataset and evolutionary model are similar to previous publications (Shimodaira and Hasegawa, 1999;Shimodaira, 2001Shimodaira, , 2002, thus allowing our proposed method to be easily compared with conventional methods. The number of unrooted trees for six taxa is 105. These trees are reordered by their likelihood values and labeled as T1, T2, . . ., T105. T1 is the ML tree as shown in Figure 2 and its tree topology is represented as (((1(23))4)56). There are three internal branches (we call them as edges) in T1, which are labeled as E1, E2, and E3. For example, E1 splits the six taxa as {23|1456} and the partition of six taxa is represented as -++---, where +/-indicates taxa 1, . . . , 6 from left to right and ++ indicates the clade {23} (we setfor taxon 6, since it is treated as the outgroup). There are 25 edges in total, and each tree is specified by selecting three edges from them, although not all the combinations of three edges are allowed.
The result of phylogenetic analysis is summarized in Table 1 for trees and Table 2 for edges. Three types of p-values are computed for each tree as well as for each edge. BP is the bootstrap probability (Felsenstein, 1985) and AU is the approximately unbiased p-value (Shimodaira, 2002). Bootstrap probabilities are computed by the non-parametric bootstrap Standard errors are shown in parentheses. Boldface indicates significance (p < 0.05) for the null hypothesis that the tree is true (outside mode). For the rest of trees (T21, . . . , T105), p-values are very small (p < 0.001). † T1 is the ML tree, i.e., the tree selected by the ML method based on the dataset of Shimodaira and Hasegawa (1999). ‡ T7 is presumably the true tree as suggested by later researches; see section 4.3.
resampling (Efron, 1979) described in section 6.1. The theory and the algorithm of BP and AU will be reviewed in section 3. Since we are testing many trees and edges at the same time, there is potentially a danger of selection bias. The issue of selection bias has been discussed in Shimodaira and Hasegawa (1999) for introducing the method of multiple comparisons of log-likelihoods (also known as Shimodaira-Hasegawa test) and in Shimodaira (2002) for introducing AU test. However, these conventional methods are only taking care of the multiplicity of comparing many log-likelihood values for computing just one p-value instead of many p-values at the same time. Therefore, we intend to further adjust the AU p-value by introducing the selective inference p-value, denoted as SI. The theory and the algorithm of SI will be explained in section 4 based on the geometric theory given in section 3. After presenting the methods, we will revisit the phyloegnetic inference in section 4.3. For developing the geometric theory in sections 3 and 4, we formulate tree selection as a mathematical formulation known as the problem of regions (Efron et al., 1996;Efron and Tibshirani, 1998). For better understanding the geometric nature of the theory, the problem of regions is explained below for phylogenetic inference, although the algorithm is simple enough to be implemented without understanding the theory. Considering the space of probability distributions (Amari and Nagaoka, 2007), the parametric models for trees are represented as manifolds in the space. The dataset (or the empirical distribution) can also be represented as a "data point" X in the space, and the ML estimates for trees are represented as projections to the manifolds. This is illustrated in the visualization of probability distributions of Figure 3A using loglikelihood vectors of models (Shimodaira, 2001), where models are simply indicated as red lines from the origin; see section 6.2 for details. This visualization may be called as model map. The point X is actually reconstructed as the minimum full model containing all the trees as submodels, and the Kullback-Leibler divergence between probability distributions is represented as the squared distance between points; see Equation (27). Computation of X is analogous to the Bayesian model averaging, but based on the ML method. For each tree, we can think of a region in the space so that this tree becomes the ML tree when X is included in the region. The regions for T1, T2, and T3 are illustrated in Figure 3B, and the region for E2 is the union of these three regions.
In Figure 3A, X is very far from any of the tree models, suggesting that all the models are wrong; the likelihood ratio statistic for testing T1 against the full model is 113.4, which is highly significant as χ 2 8 (Shimodaira, 2001, section 5). Instead of testing whether tree models are correct or not, we test whether models are significantly better than the others. As seen in Figure 3B, X is in the region for T1, meaning that the model for T1 is better than those for the other trees. For convenience, observing X in the region for T1, we state that T1 is supported by the data. Similarly, X is in the region for E2 that consists of the three regions for T1, T2, T3, thus indicating that E2 is supported by the data. Although T1 and E2 are supported by the data, there is still uncertainty as to whether the true evolutionary history of lineages is depicted because the location of X fluctuates randomly. Therefore, statistical Standard errors are shown in parentheses. Boldface without underline indicates significance (p < 0.05) for the null hypothesis that the edge is true (outside mode). Boldface with underline indicates significance (p > 0.95) for the null hypothesis that the edge is not true (inside mode). † Edges included in T1. ‡ Edges included in T7.
confidence of the outcome needs to be assessed. A mathematical procedure for statistically evaluating the outcome is provided in the following sections.

The Problem of Regions
For developing the theory, we consider (m + 1)-dimensional multivariate normal random vector Y, m ≥ 0, with unknown mean vector µ ∈ R m+1 and the identity variance matrix I m+1 : A region of interest such as tree and edge is denoted as R ⊂ R m+1 , and its complement set is denoted as There are K all regions R i , i = 1, . . . , K all , and we simply write R for one of them by dropping the index i. Observing Y = y, the null hypothesis H 0 : µ ∈ R is tested against the alternative hypothesis H 1 : µ ∈ R C . This setting is called problem of regions, and the geometric theory for non-selective inference for slightly generalized settings (e.g., exponential family of distributions) has been discussed in Efron and Tibshirani (1998) and Shimodaira (2004). This theory allows arbitrary shape of R without assuming a particular shape such as half-space or sphere, and only requires the expression (29) of section 6.3. The problem of regions is well described by geometric quantities (Figure 4). Letμ be the projection of y to the boundary surface ∂R defined asμ = arg min µ∈∂R y − µ , and β 0 be the signed distance defined as β 0 = y −μ > 0 for y ∈ R C and β 0 = − y −μ ≤ 0 for y ∈ R; see Figures 4A,B, respectively. A large β 0 indicates the evidence for rejecting H 0 : µ ∈ R, but computation of p-value will also depend on the shape of R. There should be many parameters for defining the shape, but we only need the mean curvature of ∂R atμ, which represents the amount of surface bending. It is denoted as β 1 ∈ R, and defined in (30).
Geometric quantities β 0 and β 1 of regions for trees (T1, . . . , T105) and edges (E1, . . . , E25) are plotted in Figure 5, and these values are also found in Tables 1, 2. Although the phylogenetic model of evolution for the molecular dataset X n = (x 1 , . . . , x n ) is different from the multivariate normal model (4) for y, the multiscale bootstrap method of section 3.4 estimates β 0 and β 1 using the non-parametric bootstrap probabilities (section 6.1) with bootstrap replicates X * n ′ for several values of sample size n ′ .
FIGURE 3 | Model map: Visualization of ML estimates of probability distributions for the best 15 trees. The origin represents the star-shaped tree topology (obtained by reducing the internal branches to zero length). Sites of amino acid sequences t = 1, . . . , n (black numbers) and probability distributions for trees T1, . . . ,T15 (red segments) are drawn by biplot of PCA. Auxiliary lines are drawn by hand. (A) 3-dimensional visualization using PC1, PC2, and PC3. The reconstructed data point X is also shown (green point). The ML estimates are represented as the end points of the red segments (shown by red points only for the best five trees), and they are placed on the sphere with the origin and X being placed at the poles. (B) The top-view of model map. Regions for the best three trees Ti, i = 1, 2, 3 (blue shaded regions) are illustrated; Ti will be the ML tree if X is included in the region for Ti.

Bootstrap Probability
For simulating (4) from y, we may generate replicates Y * from the bootstrap distribution ( Figure 4C) and define bootstrap probability (BP) of R as the probability of Y * being included in the region R: BP(R|y) can be interpreted as the Bayesian posterior probability P(µ ∈ R|y), because, by assuming the flat prior distribution π(µ) = constant, the posterior distribution µ|y ∼ N m+1 (y, I m+1 ) is identical to the distribution of Y * in (5). An interesting consequence of the geometric theory of Efron and Tibshirani (1998) is that BP can be expressed as where ≃ indicates the second order asymptotic accuracy, meaning that the equality is correct up to O p (n −1/2 ) with error of order O p (n −1 ); see section 6.3. For understanding the formula (7), assume that R is a half space so that ∂R is flat and β 1 = 0. Since we only have to look at the axis orthogonal to ∂R, the distribution of signed distance is identified as (1) with β 0 = z. The bootstrap distribution for (1) is Z * ∼ N(z, 1), and bootstrap probability is expressed as P(Z * ≤ 0|z) =¯ (z). Therefore, we have BP(R|y) =¯ (β 0 ). For general R with curved ∂R, the formula (7) adjusts the bias caused by β 1 . As seen in Figure 4C, R becomes smaller for β 1 > 0 than β 1 = 0, and BP becomes smaller.
BP of R C is closely related to BP of R. From the definition, The last expression also implies that the signed distance and the mean curvature of R C is −β 0 and −β 1 , respectively; this relation is also obtained by reversing the sign of v in (29).

Approximately Unbiased Test
Although BP(R|y) may work as a Bayesian confidence measure, we would like to have a frequentist confidence measure for testing H 0 : µ ∈ R against H 1 : µ ∈ R C . The signed distance of Y is denoted as β 0 (Y), and consider the region {Y | β 0 (Y) > β 0 } in which the signed distance is larger than the observed value β 0 = β 0 (y). Similar to (2), we then define an approximately unbiased (AU) p-value as where the probability is calculated for Y ∼ N m+1 (μ, I m+1 ) as illustrated in Figure 4D. The shape of the region {Y | β 0 (Y) > β 0 } is very similar to the shape of R C ; the difference is in fact only O p (n −1 ). Let us think of a point y ′ with signed distance −β 0 (shown as y in Figure 4B). Then we have where the last expression is obtained by substituting (−β 0 , β 1 ) for (β 0 , β 1 ) in (8). This formula computes AU from (β 0 , β 1 ). An intuitive interpretation of (10) is explained in section 6.4. In non-selective inference, p-values are computed using formula (10). If AU(R|y) < α, the null hypothesis H 0 : µ ∈ R is rejected and the alternative hypothesis H 1 : µ ∈ R C is accepted. This test procedure is approximately unbiased, because it controls the non-selective type-I error as and the rejection probability increases as µ moves away from R, while it decreases as µ moves into R.

Multiscale Bootstrap
In order to estimate β 0 and β 1 from bootstrap probabilities, we consider a generalization of (5) as for a variance σ 2 > 0, and define multiscale bootstrap probability of R as BP σ 2 (R|y) : = P σ 2 (Y * ∈ R|y), where P σ 2 indicates the probability with respect to (13). Although our theory is based on the multivariate normal model, the actual implementation of the algorithm uses the non-parametric bootstrap probabilities in section 6.1. To fill the gap between the two models, we consider a non-linear transformation f n so that the multivariate normal model holds at least approximately for y = f n (X n ) and Y * = f n (X * n ′ ). An example of f n is given in (25) for phylogenetic inference. Surprisingly, a specification of f n is not required for computing p-values, but we simply assume the existence of such a transformation; this property may be called as "bootstrap trick." For phylogenetic inference, we compute the non-parametric bootstrap probabilities by (24) and substitute these values for (14) with σ 2 = n/n ′ .
We can estimate β 0 and β 1 as regression coefficients by fitting the linear model (16) in terms of σ 2 to the observed values of nonparametric bootstrap probabilities (Figure 6). Interestingly, (10) is rewritten as AU(R|y) ≃¯ (ψ −1 (R|y)) by formally letting σ 2 = −1 in the last expression of (16), meaning that AU corresponds to n ′ = −n. Although σ 2 should be positive in (15), we can think of negative σ 2 in β 0 + β 1 σ 2 . See section 6.5 for details of model fitting and extrapolation to negative σ 2 .

Approximately Unbiased Test for Selective Inference
In order to argue selective inference for the problem of regions, we have to specify the selection event. Let us consider a selective region S ⊂ R m+1 so that we perform the hypothesis testing only when y ∈ S. Terada and Shimodaira (2017) considered a general shape of S, but here we treat only two special cases of S = R C and S = R; see section 6.6. Our problem is formulated as follows.
Observing Y = y from the multivariate normal model (4), we first check whether y ∈ R C or y ∈ R. If y ∈ R C and we are interested in the null hypothesis H 0 : µ ∈ R, then we may test it against the alternative hypothesis H 1 : µ ∈ R C . If y ∈ R and we are interested in the null hypothesis H 0 : µ ∈ R C , then we may test it against the alternative hypothesis H 1 : µ ∈ R. In this paper, the former case (y ∈ R C , and so β 0 > 0) is called as outside mode, and the latter case (y ∈ R, and so β 0 ≤ 0) is called as inside mode. We do not know which of the two modes of testing is performed until we observe y. Let us consider the outside mode by assuming that y ∈ R C , where β 0 > 0. Recalling that p(z, c) = p(z)/¯ (c) in section 1, we divide AU(R|y) by the selection probability to define a selective inference p-value as From the definition, SI(R|y) where BP(R C |μ) =¯ (−β 1 ) is obtained by substituting (0, β 1 ) for (β 0 , β 1 ) in (8). An intuitive justification of (18) is explained in section 6.4. For the outside mode of selective inference, p-values are computed using formula (18). If SI(R|y) < α, then reject H 0 : µ ∈ R and accept H 1 : µ ∈ R C . This test procedure is approximately unbiased, because it controls the selective type-I error as and the rejection probability increases as µ moves away from R, while it decreases as µ moves into R. Now we consider the inside mode by assuming that y ∈ R, where β 0 ≤ 0. SI of R C is obtained from (17) by exchanging the roles of R and R C .

Shortcut Computation of SI
We can compute SI from BP and AU. This will be useful for reanalyzing the results of previously published researches. Let us write BP = BP(R|y) and AU = AU(R|y). From (7) and (10), we have We can compute SI from β 0 and β 1 by (18)

Revisiting the Phylogenetic Inference
In this section, the analytical procedure outlined in section 2 is used to determine relationships among human, mouse, and rabbit. The question is: Which of mouse or human is closer to rabbit? The traditional view (Novacek, 1992) is actually supporting E6, the clade of rabbit and mouse, which is consistent with T4, T5, and T7. Based on molecular analysis, Graur et al. (1996) strongly suggested that rabbit is closer to human than mouse, thus supporting E2, which is consistent with T1, T2, and T3. However, Halanych (1998) criticized it by pointing out that E2 is an artifact caused by the long branch attraction (LBA) between mouse and opossum. In addition, Shimodaira and Hasegawa (1999) and Shimodaira (2002) suggested that T7 is not rejected by multiplicity adjusted tests. Shimodaira and Hasegawa (2005) showed that T7 becomes the ML tree by resolving the LBA using a larger dataset with more taxa. Although T1 is the ML tree based on the dataset with fewer taxa, T7 is presumably the true tree as indicated by later researches. With these observations in mind, we retrospectively interpret p-values in Tables 1, 2.
The results are shown below for the two test modes (inside and outside) as defined in section 4.1. The extent of multiplicity and selection bias depends on the number of regions under consideration, thus these numbers are considered for interpreting the results. The numbers of regions related to trees and edges are summarized in Table 3; see section 6.7 for details.
In inside mode, the null hypothesis H 0 : µ ∈ R C i is tested against the alternative hypothesis H 1 : µ ∈ R i for y ∈ R i (i.e., β 0 ≤ 0). This applies to the regions for T1, E1, E2, and E3, and they are supported by the data in the sense mentioned in the last paragraph of section 2. When H 0 is rejected by a test procedure, it is claimed that R i is significantly supported by the data, indicating H 1 holds true. For convenience, the null hypothesis H 0 is said like E1 is not true, and the alternative hypothesis H 1 is said like E1 is true; then rejection of H 0 implies that E1 is true. This procedure looks unusual, but makes sense when both R i and R C i are regions with nonzero volume. Note that selection bias can be very large in the sense that K select /K all ≈ 0 for many taxa, and non-selective tests may lead to many false positives because K true /K all ≈ 1. Therefore selective inference should be used in inside mode.
In outside mode, the null hypothesis H 0 : µ ∈ R i is tested against the alternative hypothesis H 1 : µ ∈ R C i for y ∈ R C i (i.e., β 0 > 0). This applies to the regions for T2, ..., T105, and E4, ..., E25, and they are not supported by the data. When H 0 is rejected by a test procedure, it is claimed that R i is rejected. For convenience, the null hypothesis is said like T9 is true, and the alternative hypothesis is said like T9 is not true; rejection of H 0 implies that T9 is not true. This is more or less a typical test procedure. Note that selection bias is minor in the sense that K select /K all ≈ 1 for many taxa, and non-selective tests may result in few false positives because K true /K all ≈ 0. Therefore selective inference is not much beneficial in outside mode. In addition to p-values for some trees and edges, estimated geometric quantities are also shown in the tables. We confirm that the sign of β 0 is estimated correctly for all the trees and edges. The estimated β 1 values are all positive, indicating the regions are convex. This is not surprising, because the regions are expressed as intersections of half spaces at least locally ( Figure 3B). Now p-values are examined in inside mode. (T1, E3) BP, AU, SI are all p ≤ 0.95. This indicates that T1 and E3 are not significantly supported. There are nothing claimed to be definite. (E1) BP, AU, SI are all p > 0.95, indicating E1 is significantly supported. Since E1 is associated with the best 15 trees T1, ..., T15, some of them are significantly better than the rest of trees T16, ..., T105. Significance for edges is common in phylogenetics as well as in hierarchical clustering (Suzuki and Shimodaira, 2006). (E2) The results split for this presumably wrong edge. AU > 0.95 suggests E2 is significantly supported, whereas BP, SI ≤ 0.95 are not significant. AU tends to violate the selective type-I error, leading to false positives or overconfidence in wrong trees/edges, whereas SI is approximately unbiased for the selected hypothesis. This overconfidence is explained by the inequality AU > SI (meant SI ′ here) for y ∈ R, which is obtained by comparing (12) and (20). Therefore SI is preferable to AU in inside mode. BP is safer than AU in the sense that BP < AU for β 1 > 0, but BP is not guaranteed for controlling type-I error in a frequentist sense. The two inequalities (SI, BP < AU) are verified as relative positions of the contour lines at p = 0.95 in Figure 5. The three p-values can be very different from each other for large β 1 .
Next p-values are examined in outside mode. (T2, E4, E6) BP, AU, SI are all p ≥ 0.05. They are not rejected, and there are nothing claimed to be definite. (T8, T9, ..., T105, E9,..., E25) BP, AU, SI are all p < 0.05. These trees and edges are rejected. (T7, E8) The results split for these presumably true tree and edge. BP < 0.05 suggests T7 and E8 are rejected, whereas AU, SI ≥ 0.05 are not significant. AU is approximately unbiased for controlling the type-I error when H 0 is specified in advance (Shimodaira, 2002). Since BP < AU for β 1 > 0, BP violates the type-I error, which results in overconfidence in non-rejected wrong trees. Therefore BP should be avoided in outside mode. Inequality AU < SI can be shown for y ∈ R C by comparing (10) and (18). Since the null hypothesis H 0 : µ ∈ R is chosen after looking at y ∈ R C , AU is not approximately unbiased for controlling the selective type-I error, whereas SI adjusts this selection bias. The two inequalities (BP < AU < SI) are verified as relative positions of the contour lines at p = 0.05 in Figure 5. AU and SI behave similarly (Note: K select /K all ≈ 1), while BP is very different from AU and SI for large β 1 . It is arguable which of AU and SI is appropriate: AU is preferable to SI in tree selection (K true = 1), because the multiplicity of testing is controlled as FWER = P(reject any true null) = P(AU(R true tree |Y) < α | µ ∈ R true tree ) ≤ α. The FWER is multiplied by K true ≥ 1 for edge selection, and SI does not fix it either. For testing edges in outside mode, AU may be used for screening purpose with a small α value such as α/K true .

CONCLUSION
We have developed a new method for computing selective inference p-values from multiscale bootstrap probabilities, and applied this new method to phylogenetics. It is demonstrated through theory and a real-data analysis that selective inference p-values are in particular useful for testing selected edges (i.e., clades or clusters of species) to claim that they are supported significantly if p > 1 − α. On the other hand, the previously proposed non-selective version of approximately unbiased pvalues are still useful for testing candidate trees to claim that they are rejected if p < α. Although we focused on phylogenetics, our general theory of selective inference may be applied to other model selection problems, or more general selection problems.

Bootstrap Resampling of Log-Likelihoods
Non-parametric bootstrap is often time consuming for recomputing the maximum likelihood (ML) estimates for bootstrap replicates. Kishino et al. (1990) considered the resampling of estimated log-likelihoods (RELL) method for reducing the computation. Let X n = (x 1 , . . . , x n ) be the dataset of sample size n, where x t is the site-pattern of amino acids at site t for t = 1, . . . , n. By resampling x t from X n with replacement, we obtain a bootstrap replicate X * n ′ = (x * 1 , . . . , x * n ′ ) of sample size n ′ . Although n ′ = n for the ordinary bootstrap, we will use several n ′ > 0 values for the multiscale bootstrap. The parametric model of probability distribution for tree Ti is p i (x; θ i ) for i = 1, . . . , 105, and the log-likelihood function is ℓ i (θ i ; X n ) = n t=1 log p i (x t ; θ i ). Computation of the ML estimateθ i = arg max θ i ℓ i (θ i ; X n ) is time consuming, so we do not recalculateθ * i = arg max θ i ℓ i (θ i ; X * n ′ ) for bootstrap replicates. Define the site-wise log-likelihood at site t for tree Ti as so that the log-likelihood value for tree Ti is written as ℓ i (θ i ; X n ) = n t=1 ξ ti . The bootstrap replicate of the loglikelihood value is approximated as where w * t is the number of times x t appears in X * n ′ . The accuracy of this approximation as well as the higher-order term is given in Equations (4) and (5) of Shimodaira (2001). Once ℓ i (θ * i ; X * n ′ ), i = 1, . . . , 105, are computed by (23), its ML tree is Tî * witĥ i * = arg max i=1,...,105 ℓ i (θ * i ; X * n ′ ).
The non-parametric bootstrap probability of tree Ti is obtained as follows. We generate B bootstrap replicates X * b n ′ , b = 1, . . . , B. In this paper, we used B = 10 5 . For each X * b n ′ , the ML tree Tî * b is computed by the method described above. Then we count the frequency that Ti becomes the ML tree in the B replicates. The non-parametric bootstrap probability of tree Ti is computed by The non-parametric bootstrap probability of a edge is computed by summing BP(Ti, n ′ ) over the associated trees. An example of the transformation Y * = f n (X * n ′ ) mentioned in section 3.4 is where L * and V n is the variance matrix of L * n . According to the approximation (23) and the central limit theorem, (13) holds well for sufficiently large n and n ′ with m = 104 and σ 2 = n/n ′ . It also follows from the above argument that var(ℓ * i − ℓ * j ) ≈ (n ′ /n) ξ i − ξ j 2 , and thus the variance of log-likelihood difference is which gives another insight into the visualization of section 6.2, where the variance can be interpreted as the divergence between the two models; see Equation (27). This approximation holds well when the two predictive distributions p i (x;θ i ), p j (x;θ j ) are not very close to each other. When they are close to each other, however, the higher-order term ignored in (26) becomes dominant, and there is a difficulty for deriving the limiting distribution of the log-likelihood difference in the model selection test (Shimodaira, 1997;Schennach and Wilhelm, 2017).

Visualization of Probability Models
For representing the probability distribution of tree Ti, we define ξ i : = (ξ 1i , . . . , ξ ni ) T ∈ R n from (22) for i = 1, . . . , 15. The idea behind the visualization of Figure 3 is that locations of ξ i in R n will represent locations of p i (x;θ i ) in the space of probability distributions. Let D KL (p i p j ) be the Kullback-Leibler divergence between the two distributions. For sufficiently small (1/n) ξ i − ξ j 2 , the squared distance in R n approximates n times Jeffreys divergence for non-nested models (Shimodaira, 2001, section 6). When a model p 0 is nested in p i , it becomes ξ i − ξ 0 2 ≈ 2n × D KL (p i (x;θ i ) p 0 (x;θ 0 )) ≈ 2 × (ℓ i (θ i ; X n ) − ℓ 0 (θ 0 ; X n )). We explain three different visualizations of Figure 7. There are only minor differences between the plots, and the visualization is not sensitive to the details. For dimensionality reduction, we have to specify the origin c ∈ R n and consider vectors a i : = ξ i −c. A naive choice would be the average c = 15 i=1 ξ i /15. By applying PCA without centering and scaling (e.g., prcomp with option center=FALSE, scale=FALSE in R) to the matrix (a 1 , . . . , a 15 ), we obtain the visualization of ξ i as the axes (red arrows) of biplot in Figure 7A.
For computing the "data point" X in Figure 3, we need more models. Let tree T106 be the star topology with no internal branch (completely unresolved tree), and T107, . . . , T131 be partially resolved tree topologies with only one internal branch corresponding to E1, . . . , E25, whereas T1, . . . , T105 are fully resolved trees (bifurcating trees). Then define η i : = ξ 106+i , i = 0, . . . , 25. Now we take c = η 0 for computing a i = ξ i − η 0 and b i = η i − η 0 . There is hierarchy of models: η 0 is the submodel nested in all the other models, and η 1 , η 2 , η 3 , for example, are submodels of ξ 1 (T1 includes E1, E2, E3). By combining these non-nested models, we can reconstruct a comprehensive model in which all the other models are nested as submodels (Shimodaira, 2001, Equation 10 in section 5). The idea is analogous to reconstructing the full model y = β 1 x 1 + · · · + β 25 x 25 + ǫ of multiple regression from submodels y = β 1 x 1 + ǫ, . . . , y = β 25 x 25 + ǫ. Thus we call it as "full model" in this paper, and the ML estimate of the full model is indicated as the data point X; it is also said "super model" in Shimodaira and Hasegawa (2005). Let B = (b 1 , . . . , b 25 ) ∈ R n×25 and d = ( b 1 2 , . . . , b 25 2 ) T ∈ R 25 , then the vector for the full model is computed approximately by For the visualization of the best 15 trees, we may use only b 1 , . . . , b 11 , because they include E1 and two more edges from E2, . . . ,E11. In Figures 3, 7B, we actually modified the above computation slightly so that the star topology T106 is replaced by T107, the partially resolved tree corresponding to E1 (T107 is also said star topology by treating clade (23) as a leaf of the tree), and the 10 partially resolved trees for E2, . . . , E11 are replaced by those for (E1,E2), . . . , (E1,E11), respectively; the origin becomes the maximal model nested in all the 15 trees, and X becomes the minimal full model containing all the 15 trees. Just before applying PCA in Figure 7B, a 1 , . . . , a 15 are projected to the space orthogonal to a X , so that the plot becomes the "top-view" of Figure 3A with a X being at the origin.
In Figure 7C, we attempted a even simpler computation without using ML estimates for partially resolved trees. We used B = (a 1 , . . . , a 15 ) and d = ( a 1 2 , . . . , a 15 2 ) T , and taking the largest 10 singular values for computing the inverse in (28). The orthogonal projection to a X is applied before PCA.

Asymptotic Theory of Smooth Surfaces
For expressing the shape of the region R ⊂ R m+1 , we use a local coordinate system (u, v) ∈ R m+1 with u ∈ R m , v ∈ R. In a neighborhood of y, the region is expressed as so that y = (0, β 0 ) (i.e., u = (0, . . . , 0) and v = β 0 ), and h(0) = 0, ∂h/∂u i | 0 = 0, i = 1, . . . , m. The projection now becomes the originμ = (0, 0), and the signed distance is β 0 . The mean curvature of surface ∂R atμ is now defined as which is interpreted as the trace of the hessian matrix of h. When R is convex at least locally in the neighborhood, all the eigenvalues of the hessian are non-negative, leading to β 1 ≥ 0, whereas concave R leads to β 1 ≤ 0. In particular, β 1 = 0 when ∂R is flat (i.e., h(u) ≡ 0). Since the transformation y = f n (X n ) depends on n, the shape of the region R actually depends on n, although the dependency is implicit in the notation. As n goes larger, the standard deviation of estimates, in general, reduces at the rate n −1/2 . For keeping the variance constant in (4), we actually magnifying the space by the factor n 1/2 , meaning that the boundary surface ∂R approaches flat as n → ∞. More specifically, the magnitude of mean curvature is of order β 1 = O p (n −1/2 ). The magnitude of ∂ 3 h/∂u i ∂u j ∂u k and higher order derivatives is O p (n −1 ), and we ignore these terms in our asymptotic theory. For keeping µ = O(1) in (4), we also consider the setting of "local alternatives, " meaning that the parameter values approach a origin on the boundary at the rate n −1/2 .

Bridging the Problem of Regions to the Z-Test
Here we explain the problem of regions in terms of the ztest by bridging the multivariate problem of section 3 to the 1-dimensional case of section 1.

Model Fitting in Multiscale Bootstrap
We have used thirteen σ 2 values from 1/9 to 9 (equally spaced in log-scale). This range is relatively large, and we observe a slight deviation from the linear model β 0 +β 1 σ 2 in Figure 6. Therefore we fit other models to the observed values of ψ σ 2 as implemented in scaleboot package (Shimodaira, 2008). For example, poly.k model is k−1 i=0 β i σ 2i , and sing.3 model is β 0 + β 1 σ 2 (1 + β 2 (σ − 1)) −1 . In Figure 6A, poly.3 is the best model according to AIC (Akaike, 1974). In Figure 6B, poly.2, poly.3, and sing.3 are combined by model averaging with Akaike weights. Then β 0 and β 1 are estimated from the tangent line to the fitted curve of ψ σ 2 at σ 2 = 1. In Figure 6, the tangent line is drawn as red line for extrapolating ψ σ 2 to σ 2 = −1. Shimodaira (2008) and Terada and Shimodaira (2017) considered the Taylor expansion of ψ σ 2 at σ 2 = 1 as a generalization of the tangent line for improving the accuracy of AU and SI.
In the implementation of CONSEL (Shimodaira and Hasegawa, 2001) and pvclust (Suzuki and Shimodaira, 2006), we use a narrower range of σ 2 values (ten σ −2 values: 0.5, 0.6, . . . , 1.4). Only the linear model β 0 + β 1 σ 2 is fitted there. The estimated β 0 and β 1 should be very close to those estimated from the tangent line described above. An advantage of using wider range of σ 2 in scaleboot is that the standard error of β 0 and β 1 will become smaller.

General Formula of Selective Inference
Let H, S ⊂ R m+1 be regions for the null hypothesis and the selection event, respectively. We would like to test the null hypothesis H 0 : µ ∈ H against the alternative H 1 : µ ∈ H C conditioned on the selection event y ∈ S. We have considered the outside mode H = R, S = R C in (18) and the inside mode (20). For a general case of H, S, Terada and Shimodaira (2017) gave a formula of approximately unbiased p-value of selective inference as where geometric quantities β 0 , β 1 are defined for the regions H, S. We assumed that H and S C are expressed as (29), and two surfaces ∂H, ∂S are nearly parallel to each other with tangent planes differing only O p (n −1/2 ). The last assumption always holds for (18), because ∂H = ∂R and ∂S = ∂R C are identical and of course parallel to each other.
Here we explain why we have considered the special case of S = H C for phylogenetic inference. First, we suppose that the selection event satisfies S ⊂ H C , because a reasonable test would not reject H 0 unless y ∈ H C . Note that y ∈ S ⊂ H C implies 0 ≤ −β S 0 ≤ β H 0 . Therefore, β H 0 + β S 0 ≥ 0 leads to SI(H|S, y) ≥ SI(H|y), where SI(H|y) : = SI(H|H C , y) is obtained from (33) by letting β H 0 + β S 0 = 0 for S = H C . The p-value SI(H|S, y) becomes smaller as S grows, and S = H C gives the smallest p-value, leading to the most powerful selective test. Therefore the choice S = H C is preferable to any other choice of selection event satisfying S ⊂ H C . This kind of property is mentioned in Fithian et al. (2014) as the monotonicity of selective error in the context of "data curving." Let us see how these two p-values differ for the case of E2 by specifying H = R C E2 and S = R T1 . In this case, the two surfaces ∂H, ∂S may not be very parallel to each other, thus violating the assumption of SI(H|S, y), so we only intend to show the potential difference between the two p-values. The geometric quantities are where SI ′ (R C E2 |y) = 1 − SI(R C E2 |y) = 0.903 is shown in Table 2. As you see, SI(H|y) is easier to reject H 0 than SI(H|S, y).

Number of Regions for Phylogenetic Inference
The regions R i , i = 1, . . . , K all correspond to trees or edges. In inside and outside modes, the number of total regions is K all = 105 for trees and K all = 25 for edges when the number of taxa is N = 6. For general N ≥ 3, they grow rapidly as K all = (2N − 5)!/(2 N−3 (N − 3)!) for trees and K all = 2 N−1 − (N + 1) for edges. Next consider the number of selected regions K select . In inside mode, regions with y ∈ R i are selected, and the number is counted as K select = 1 for trees and K select = N − 3 = 3 for edges. In outside mode, regions with y ∈ R i are selected, and thus the number is K all minus that for inside mode; K select = K all − 1 = 104 for trees and K select = K all − (N − 3) = 22 for edges. Finally, consider the number of true null hypotheses, denoted as K true . The null hypothesis holds true when µ ∈ R i in inside mode and µ ∈ R i in outside mode, and thus K true is the same as the number of regions with y ∈ R i in inside mode and y ∈ R i in outside mode (These numbers do not depend on the value of y by ignoring the case of y ∈ ∂R i ). Therefore, K true = K all − K select for both cases.

Selective Inference of Lasso Regression
Selective inference is considered for the variable selection of regression analysis. Here, we deal with prostate cancer data (Stamey et al., 1989) in which we predict the level of prostatespecific antigen (PSA) from clinical measures. The dataset is available in the R package ElemStatLearn (Halvorsen, 2015). We consider a linear model to the log of PSA (lpsa), with 8 predictors such as the log prostate weight (lweight), age, and so on. All the variables are standardized to have zero mean and unit variance.
The goal is to provide the valid selective inference for the partial regression coefficients of the selected variables by lasso (Tibshirani, 1996). Let n and p be the number of observations and the number of predictors.M is the set of selected variables, andŝ represents the signs of the selected regression coefficients. We suppose that regression responses are distributed as Y ∼ N(µ, τ 2 I n ) where µ ∈ R n and τ > 0. Let e i be the ith residual. Resampling the scaled residuals σ e i (i = 1, . . . , n) with several values of scale σ 2 , we can apply the multiscale bootstrap method described in section 4 for the selective inference in the regression problem. Here, we note that the target of the inference is the true partial regression coefficients: where X ∈ R n×p is the design matrix. We compute four types of intervals with confidence level 1 − α = 0.95 for selected variable j. is the selective confidence interval under the selected model proposed by Lee et al. (2016) and Tibshirani et al. (2016), which is computed by fixedLassoInf with type="full" in R package selectiveInference (Tibshirani et al., 2017) ], which is the selective confidence interval under the selection event that variable j is selected. These three confidence intervals are exact, in the sense that Note that the selection event of variable j, i.e., {j ∈M,ŝ j } can be represented as a union of polyhedra on R n , and thus, according to the polyhedral lemma (Lee et al., 2016;Tibshirani et al., 2016), we can compute a valid confidence interval [L variable j , U variable j ]. However, this computation is prohibitive for p > 10, because all the possible combinations of models with variable j are considered. Therefore, we compute its approximation [L variable j ,Û variable j ] by the multiscale bootstrap method of section 4 with much faster computation even for larger p.
The confidence intervals are shown in Figure 1. For adjusting the selection bias, the three confidence intervals of selective inference are longer than the ordinary confidence interval. Comparing [L model j , U model j ] and [L variable j , U variable j ], the latter is shorter, and would be preferable. This is because the selection event of the latter is less restrictive as {M,ŝ} ⊆ {j ∈M,ŝ j }; see section 6.6 for the reason why larger selection event is better.

AUTHOR CONTRIBUTIONS
HS and YT developed the theory of selective inference. HS programmed the multiscale bootstrap software and conducted the phylogenetic analysis. YT conducted the lasso analysis. HS wrote the manuscript. All authors have approved the final version of the manuscript.

FUNDING
This research was supported in part by JSPS KAKENHI Grant (16H02789 to HS, 16K16024 to YT).