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−α, whereas the non-selective p-value is still useful for testing candidate trees to claim that they are rejected if p < α. 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 model of scaling-law to the non-parametric multiscale bootstrap probabilities. 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.

1. 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 θ ∈ ℝ and variance 1:

$\begin{array}{ll}Z~N(\theta ,1).& (1)\end{array}$

Observing Z = z, we would like to test the null hypothesis H₀ : θ ≤ 0 against the alternative hypothesis H₁ : θ > 0. We denote the cumulative distribution function of N(0, 1) as Φ(x) and define the upper tail probability as $\bar{\Phi }(x)=1-\Phi (x)=\Phi (-x)$. Then, the ordinary (i.e., non-selective) inference leads to the p-value of the one-tailed z-test as

$\begin{array}{ll}p(z):=P(Z>z\mid \theta =0)=\bar{\Phi }(z).& (2)\end{array}$

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)

$\begin{array}{ll}p(z,c):=P(Z>z\mid Z>c,\theta =0)=\bar{\Phi }(z)/\bar{\Phi }(c),& (3)\end{array}$

where p(z) of Equation (2) is divided by the selection probability $P(Z>c\mid \theta =0)=\bar{\Phi }(c)$. In the case of c = 0, this corresponds to the two-tailed z-test, because the selection probability is $\bar{\Phi }(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)<\alpha \mid Z>c,\theta =0)=\alpha /\bar{\Phi }(c)>\alpha$. The selection bias can be very large when $\bar{\Phi }(c)\ll 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, 2004, 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.

FIGURE 1

Figure 1. Confidence intervals of regression coefficients for selected variables by lasso; see section 6.8 for details. All intervals are computed for confidence level 1 − α at α = 0.01. (Black) the ordinary confidence interval $[L_j^{\mathrm{\text{ordinary}}},U_j^{\mathrm{\text{ordinary}}}]$. (Green) the selective confidence interval $[L_j^{\mathrm{\text{model}}},U_j^{\mathrm{\text{model}}}]$ under the selected model. (Blue) the selective confidence interval $[L_j^{\mathrm{\text{variable}}},U_j^{\mathrm{\text{variable}}}]$ under the selection event that variable j is selected. (Red) the multiscale bootstrap version of selective confidence interval $[{\widehat{L}}_j^{\mathrm{\text{variable}}},{\widehat{U}}_j^{\mathrm{\text{variable}}}]$ under the selection event that variable j is selected.

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 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.

2. 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,\dots ,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, 2001, 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 set - for 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.

FIGURE 2

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.

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 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.

TABLE 1

Table 1. Three types of p-values (BP, AU, SI) and geometric quantities (β₀, β₁) for the best 20 trees.

TABLE 2

Table 2. Three types of p-values (BP, AU, SI) and geometric quantities (β₀, β₁) for all the 25 edges of six taxa.

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 log-likelihood 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.

FIGURE 3

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.

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 ${\chi }_8^2$ (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 confidence of the outcome needs to be assessed. A mathematical procedure for statistically evaluating the outcome is provided in the following sections.

3. Non-Selective Inference for the Problem of Regions

3.1. The Problem of Regions

For developing the theory, we consider (m + 1)-dimensional multivariate normal random vector Y, m ≥ 0, with unknown mean vector μ ∈ ℝ^m+1 and the identity variance matrix I_m+1:

$\begin{array}{ll}Y~N_{m+1}(\mu ,I_{m+1}).& (4)\end{array}$

A region of interest such as tree and edge is denoted as ${R}\subset {ℝ}^{m+1}$, and its complement set is denoted as ${{R}}^C={ℝ}^{m+1}\backslash {R}$. 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:\mu \in {R}$ is tested against the alternative hypothesis $H_1:\mu \in {{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 $\widehat{\mu }$ be the projection of y to the boundary surface $\partial {R}$ defined as

$\begin{array}{l}\widehat{\mu }={\displaystyle \underset{\mu \in \partial {R}}{arg min}}\Vert y-\mu \Vert ,\end{array}$

and β₀ be the signed distance defined as ${\beta }_0=\Vert y-\widehat{\mu }\Vert >0$ for $y\in {{R}}^C$ and ${\beta }_0=-\Vert y-\widehat{\mu }\Vert \le 0$ for $y\in {R}$; see Figures 4A,B, respectively. A large β₀ indicates the evidence for rejecting $H_0:\mu \in {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 $\partial {R}$ at $\widehat{\mu }$, which represents the amount of surface bending. It is denoted as β₁ ∈ ℝ, and defined in (30).

FIGURE 4

Figure 4. Problem of regions. (A) β₀ > 0 when $y\in {{R}}^C$, then select the null hypothesis $\mu \in {R}$. (B) β₀ ≤ 0 when $y\in {R}$, then select the null hypothesis $\mu \in {{R}}^C$. (C) The bootstrap distribution of $y^{*}~N_{m+1}(y,I_{m+1})$ (red shaded distribution). (D) The null distribution of $Y~N_{m+1}(\widehat{\mu },I_{m+1})$ (green shaded distribution).

Geometric quantities β₀ and β₁ 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,\dots ,x_n)$ is different from the multivariate normal model (4) for y, the multiscale bootstrap method of section 3.4 estimates β₀ and β₁ using the non-parametric bootstrap probabilities (section 6.1) with bootstrap replicates ${{X}}_{n^{\prime }}^{*}$ for several values of sample size n′.

FIGURE 5

Figure 5. Geometric quantities of regions (β₀ and β₁) for trees and edges are estimated by the multiscale bootstrap method (section 3.4). The three types of p-value (BP, AU, SI) are computed from β₀ and β₁, and their contour lines are drawn at p = 0.05 and 0.95.

3.2. Bootstrap Probability

For simulating (4) from y, we may generate replicates Y^* from the bootstrap distribution (Figure 4C)

$\begin{array}{ll}{Y}^{*}~N_{m+1}(y,I_{m+1}),& (5)\end{array}$

and define bootstrap probability (BP) of ${R}$ as the probability of Y^* being included in the region ${R}$:

$\begin{array}{ll}\text{BP}({R}|y):=P(Y^{*}\in {R}|y).& (6)\end{array}$

$\mathrm{\text{BP}}({R}|y)$ can be interpreted as the Bayesian posterior probability $P(\mu \in {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

$\begin{array}{ll}\text{BP}({R}|y)\simeq \bar{\Phi }({\beta }_0+{\beta }_1),& (7)\end{array}$

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 $\partial {R}$ is flat and β₁ = 0. Since we only have to look at the axis orthogonal to $\partial {R}$, the distribution of signed distance is identified as (1) with β₀ = z. The bootstrap distribution for (1) is Z^* ~ N(z, 1), and bootstrap probability is expressed as $P(Z^{*}\le 0|z)=\bar{\Phi }(z)$. Therefore, we have $\mathrm{\text{BP}}({R}|y)=\bar{\Phi }({\beta }_0)$. For general ${R}$ with curved $\partial {R}$, the formula (7) adjusts the bias caused by β₁. As seen in Figure 4C, ${R}$ becomes smaller for β₁ > 0 than β₁ = 0, and BP becomes smaller.

BP of ${{R}}^C$ is closely related to BP of ${R}$. From the definition,

$\begin{array}{ll}\text{BP}({{R}}^C|y)=1-\text{BP}({R}|y)\simeq 1-\bar{\Phi }({\beta }_0+{\beta }_1)=\bar{\Phi }(-{\beta }_0-{\beta }_1).& (8)\end{array}$

The last expression also implies that the signed distance and the mean curvature of ${{R}}^C$ is −β₀ and −β₁, respectively; this relation is also obtained by reversing the sign of v in (29).

3.3. Approximately Unbiased Test

Although $\mathrm{\text{BP}}({R}|y)$ may work as a Bayesian confidence measure, we would like to have a frequentist confidence measure for testing $H_0:\mu \in {R}$ against $H_1:\mu \in {{R}}^C$. The signed distance of Y is denoted as β₀(Y), and consider the region {Y ∣ β₀(Y) > β₀} in which the signed distance is larger than the observed value β₀ = β₀(y). Similar to (2), we then define an approximately unbiased (AU) p-value as

$\begin{array}{ll}\text{AU}({R}|y):=P({\beta }_0(Y)>{\beta }_0\mid \mu =\widehat{\mu })=\text{BP}(\{Y\mid {\beta }_0(Y)>{\beta }_0\}|\widehat{\mu }),& (9)\end{array}$

where the probability is calculated for $Y~N_{m+1}(\widehat{\mu },I_{m+1})$ as illustrated in Figure 4D. The shape of the region {Y ∣ β₀(Y) > β₀} 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 −β₀ (shown as y in Figure 4B). Then we have

$\begin{array}{ll}\text{AU}({R}|y)\simeq \text{BP}({{R}}^C|y^{\prime })\simeq \bar{\Phi }({\beta }_0-{\beta }_1),& (10)\end{array}$

where the last expression is obtained by substituting (−β₀, β₁) for (β₀, β₁) in (8). This formula computes AU from (β₀, β₁). An intuitive interpretation of (10) is explained in section 6.4.

In non-selective inference, p-values are computed using formula (10). If $\mathrm{\text{AU}}({R}|y)<\alpha$, the null hypothesis $H_0:\mu \in {R}$ is rejected and the alternative hypothesis $H_1:\mu \in {{R}}^C$ is accepted. This test procedure is approximately unbiased, because it controls the non-selective type-I error as

$\begin{array}{ll}P\left(\text{AU}({R}|y)<\alpha \mid \mu \in \partial {R}\right)\simeq \alpha ,& (11)\end{array}$

and the rejection probability increases as μ moves away from ${R}$, while it decreases as μ moves into ${R}$.

Exchanging the roles of ${R}$ and ${{R}}^C$ also allows for another hypothesis testing. AU of ${{R}}^C$ is obtained from (9) by reversing the inequality as $\mathrm{\text{AU}}({{R}}^C|y)=\mathrm{\text{BP}}(\{Y\mid {\beta }_0(Y)<{\beta }_0\}|\widehat{\mu })=1-\mathrm{\text{AU}}({R}|y)$. This is also confirmed by substituting (−β₀, −β₁), i.e., the geometric quantities of ${{R}}^C$, for (β₀, β₁) in (10) as

$\begin{array}{ll}\text{AU}({{R}}^C|y)\simeq \bar{\Phi }(-{\beta }_0+{\beta }_1)\simeq 1-\text{AU}({R}|y).& (12)\end{array}$

If $\mathrm{\text{AU}}({{R}}^C|y)<\alpha$ or equivalently $\mathrm{\text{AU}}({R}|y)>1-\alpha$, then we reject $H_0:\mu \in {{R}}^C$ and accept $H_1:\mu \in {R}$.

3.4. Multiscale Bootstrap

In order to estimate β₀ and β₁ from bootstrap probabilities, we consider a generalization of (5) as

$\begin{array}{ll}{Y}^{*}~N_{m+1}(y,{\sigma }^2{I}_{m+1}),& (13)\end{array}$

for a variance σ² > 0, and define multiscale bootstrap probability of ${R}$ as

$\begin{array}{ll}{\text{BP}}_{{\sigma }^2}({R}|y):=P_{{\sigma }^2}(Y^{*}\in {R}|y),& (14)\end{array}$

where $P_{{\sigma }^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^{\prime }}^{*})$. 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 σ² = n/n′.

For estimating β₀ and β₁, we need to have a scaling law which explains how ${\mathrm{\text{BP}}}_{{\sigma }^2}$ depends on the scale σ. We rescale (13) by multiplying σ⁻¹ so that ${\sigma }^{-1}{Y}^{*}~N_{m+1}({\sigma }^{-1}y,I_{m+1})$ has the variance σ² = 1. y and ${R}$ are now resaled by the factor σ⁻¹, which amounts to signed distance ${\beta }_0{\sigma }^{-1}$ and mean curvature β₁σ (Shimodaira, 2004). Therefore, by substituting $({\beta }_0{\sigma }^{-1},{\beta }_1\sigma )$ for (β₀, β₁) in (7), we obtain

$\begin{array}{ll}{\text{BP}}_{{\sigma }^2}({R}|y)\simeq \bar{\Phi }({\beta }_0{\sigma }^{-1}+{\beta }_1\sigma ).& (15)\end{array}$

For better illustrating how ${\mathrm{\text{BP}}}_{{\sigma }^2}$ depends on σ², we define

$\begin{array}{ll}{\psi }_{{\sigma }^2}({R}|y):=\sigma {\bar{\Phi }}^{-1}({\text{BP}}_{{\sigma }^2}({R}|y))\simeq {\beta }_0+{\beta }_1{\sigma }^2.& (16)\end{array}$

We can estimate β₀ and β₁ as regression coefficients by fitting the linear model (16) in terms of σ² to the observed values of non-parametric bootstrap probabilities (Figure 6). Interestingly, (10) is rewritten as $\mathrm{\text{AU}}({R}|y)\simeq \bar{\Phi }({\psi }_{-1}({R}|y))$ by formally letting σ² = −1 in the last expression of (16), meaning that AU corresponds to n′ = −n. Although σ² should be positive in (15), we can think of negative σ² in ${\beta }_0+{\beta }_1{\sigma }^2$. See section 6.5 for details of model fitting and extrapolation to negative σ².

FIGURE 6

Figure 6. Multiscale bootstrap for (A) tree T1 and (B) edge E2. ${\psi }_{{\sigma }^2}({R}|y)$ is computed by the non-parametric bootstrap probabilities for several σ² = n/n′ values, then β₀ and β₁ are estimated as the intercept and the slope, respectively. See section 6.5 for details.

4. Selective Inference for the Problem of Regions

4.1. 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}\subset {{R}}^{m+1}$ so that we perform the hypothesis testing only when $y\in {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\in {{R}}^C$ or $y\in {R}$. If $y\in {{R}}^C$ and we are interested in the null hypothesis $H_0:\mu \in {R}$, then we may test it against the alternative hypothesis $H_1:\mu \in {{R}}^C$. If $y\in {R}$ and we are interested in the null hypothesis $H_0:\mu \in {{R}}^C$, then we may test it against the alternative hypothesis $H_1:\mu \in {R}$. In this paper, the former case ($y\in {{R}}^C$, and so β₀ > 0) is called as outside mode, and the latter case ($y\in {R}$, and so β₀ ≤ 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\in {{R}}^C$, where β₀ > 0. Recalling that $p(z,c)=p(z)/\bar{\Phi }(c)$ in section 1, we divide $\mathrm{\text{AU}}({R}|y)$ by the selection probability to define a selective inference p-value as

$\begin{array}{ll}\text{SI}({R}|y):=\frac{P({\beta }_0(Y)>{\beta }_0\mid \mu =\widehat{\mu })}{P(Y\in {{R}}^C\mid \mu =\widehat{\mu })}=\frac{\text{AU}({R}|y)}{\text{BP}({{R}}^C|\widehat{\mu })}.& (17)\end{array}$

From the definition, $\mathrm{\text{SI}}({R}|y)\in (0,1)$, because $\left\{Y\mid {\beta }_0(Y)>{\beta }_0\right\}\subset {{R}}^C$ for β₀ > 0. This p-value is computed from (β₀, β₁) by

$\begin{array}{ll}\text{SI}({R}|y)\simeq \frac{\bar{\Phi }({\beta }_0-{\beta }_1)}{\bar{\Phi }(-{\beta }_1)},& (18)\end{array}$

where $\mathrm{\text{BP}}({{R}}^C|\widehat{\mu })=\bar{\Phi }(-{\beta }_1)$ is obtained by substituting (0, β₁) for (β₀, β₁) 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 $\mathrm{\text{SI}}({R}|y)<\alpha$, then reject $H_0:\mu \in {R}$ and accept $H_1:\mu \in {{R}}^C$. This test procedure is approximately unbiased, because it controls the selective type-I error as

$\begin{array}{ll}P\left(\text{SI}({R}|Y)<\alpha \mid Y\in {{R}}^C,\mu \in \partial {R}\right)\simeq \alpha ,& (19)\end{array}$

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\in {R}$, where β₀ ≤ 0. SI of ${{R}}^C$ is obtained from (17) by exchanging the roles of ${R}$ and ${{R}}^C$.

$\begin{array}{ll}\text{SI}\left({{R}}^C|y\right)=\frac{\text{AU}\left({{R}}^C|y\right)}{\text{BP}\left({R}|\widehat{\mu }\right)}\simeq \frac{\overline{\Phi }\left(-{\beta }_0+{\beta }_1\right)}{\overline{\Phi }\left({\beta }_1\right)}.& (20)\end{array}$

For the inside mode of selective inference, p-values are computed using formula (20). If $\mathrm{\text{SI}}({{R}}^C|y)<\alpha$, then reject $H_0:\mu \in {{R}}^C$ and accept $H_1:\mu \in {R}$. Unlike the non-selective p-value $\mathrm{\text{AU}}({{R}}^C|y)$, $\mathrm{\text{SI}}({{R}}^C|y)<\alpha$ is not equivalent to $\mathrm{\text{SI}}({R}|y)>1-\alpha$, because $\mathrm{\text{SI}}({R}|y)+\mathrm{\text{SI}}({{R}}^C|y)\ne 1$. For convenience, we define

$\begin{array}{ll}{\text{SI}}^{\prime }({R}|y):=\{\begin{array}{cc}\text{SI}({R}|y)\hfill & y\in {{R}}^C\hfill \\ 1-\text{SI}({{R}}^C|y)\hfill & y\in {R}\hfill \end{array}\hfill & (21)\end{array}$

so that SI′ > 1 − α implies $\mathrm{\text{SI}}({{R}}^C|y)<\alpha$. In our numerical examples of Figure 5, Tables 1, 2, SI′ is simply denoted as SI. We do not need to consider (21) for BP and AU, because ${\mathrm{\text{BP}}}^{\prime }({R}|y)=\mathrm{\text{BP}}({R}|y)$ and ${\mathrm{\text{AU}}}^{\prime }({R}|y)=\mathrm{\text{AU}}({R}|y)$ from (8) and (12).

4.2. 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 $\mathrm{\text{BP}}=\mathrm{\text{BP}}({R}|y)$ and $\mathrm{\text{AU}}=\mathrm{\text{AU}}({R}|y)$. From (7) and (10), we have

$\begin{array}{l}{\beta }_0=\frac{1}{2}\left({\bar{\Phi }}^{-1}(\text{BP})+{\bar{\Phi }}^{-1}(\text{AU})\right)\\ {\beta }_1=\frac{1}{2}\left({\bar{\Phi }}^{-1}(\text{BP})-{\bar{\Phi }}^{-1}(\text{AU})\right).\end{array}$

We can compute SI from β₀ and β₁ by (18) or (20). More directly, we may compute

$\begin{array}{l}\text{SI}({R}|y)=\frac{\text{AU}}{\bar{\Phi }\left\{\frac{1}{2}\left({\bar{\Phi }}^{-1}(\text{AU})-{\bar{\Phi }}^{-1}(\text{BP})\right)\right\}}\\ \text{SI}({{R}}^C|y)=\frac{1-\text{AU}}{\bar{\Phi }\left\{\frac{1}{2}\left({\bar{\Phi }}^{-1}(\text{BP})-{\bar{\Phi }}^{-1}(\text{AU})\right)\right\}}.\end{array}$

4.3. 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.

TABLE 3

Table 3. The number of regions for trees and edges. The number of taxa is N = 6.

In inside mode, the null hypothesis $H_0:\mathit{\mu }\in {{R}}_i^C$ is tested against the alternative hypothesis $H_1:\mathit{\mu }\in {{R}}_i$ for $y\in {{R}}_i$ (i.e., β₀ ≤ 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₀ is rejected by a test procedure, it is claimed that ${{R}}_i$ is significantly supported by the data, indicating H₁ holds true. For convenience, the null hypothesis H₀ is said like E1 is not true, and the alternative hypothesis H₁ is said like E1 is true; then rejection of H₀ implies that E1 is true. This procedure looks unusual, but makes sense when both ${{R}}_i$ and ${{R}}_i^C$ 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:\mathit{\mu }\in {{R}}_i$ is tested against the alternative hypothesis $H_1:\mathit{\mu }\in {{R}}_i^C$ for $y\in {{R}}_i^C$ (i.e., β₀ > 0). This applies to the regions for T2, …, T105, and E4, …, E25, and they are not supported by the data. When H₀ 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₀ 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 β₀ is estimated correctly for all the trees and edges. The estimated β₁ 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\in {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 β₁ > 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 β₁.

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₀ is specified in advance (Shimodaira, 2002). Since BP < AU for β₁ > 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\in {{R}}^C$ by comparing (10) and (18). Since the null hypothesis $H_0:\mathit{\mu }\in {R}$ is chosen after looking at $y\in {{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 β₁. 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 $\mathrm{\text{FWER}}=P(\mathrm{\text{reject any true null}})=P(\mathrm{\text{AU}}({{R}}_{\mathrm{\text{true tree}}}|y)<\alpha \mid \mathit{\mu }\in {{R}}_{\mathrm{\text{true tree}}})\le \alpha$. 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.

5. 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 p-values 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.

6. Remarks

6.1. 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,\dots ,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^{\prime }}^{*}=(x_1^{*},\dots ,x_{n^{\prime }}^{*})$ 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 ${\ell }_i({\theta }_i;{{X}}_n)=\sum _{t=1}^{n}log{p}_i(x_t;{\theta }_i)$. Computation of the ML estimate ${\widehat{\theta }}_i={arg max}_{{\theta }_i}{\ell }_i({\theta }_i;{{X}}_n)$ is time consuming, so we do not recalculate ${\widehat{\theta }}_i^{*}={arg max}_{{\theta }_i}{\ell }_i({\theta }_i;{{X}}_{n^{\prime }}^{*})$ for bootstrap replicates. Define the site-wise log-likelihood at site t for tree Ti as

$\begin{array}{ll}{\xi }_{ti}=log{p}_i(x_t;{\widehat{\theta }}_i),t=1,\dots ,n,i=1,\dots ,105,& (22)\end{array}$

so that the log-likelihood value for tree Ti is written as ${\ell }_i({\widehat{\theta }}_i;{{X}}_n)=\sum _{t=1}^n{\xi }_{ti}$. The bootstrap replicate of the log-likelihood value is approximated as

$\begin{array}{ll}{\ell }_i({\widehat{\theta }}_i^{*};{{X}}_{n^{\prime }}^{*})\approx {\ell }_i({\widehat{\theta }}_i;{{X}}_{n^{\prime }}^{*})={\displaystyle \sum _{t=1}^n}{w}_t^{*}{\xi }_{ti},& (23)\end{array}$

where $w_t^{*}$ is the number of times x_t appears in ${{X}}_{n^{\prime }}^{*}$. The accuracy of this approximation as well as the higher-order term is given in Equations (4) and (5) of Shimodaira (2001). Once ${\ell }_i({\widehat{\theta }}_i^{*};{{X}}_{n^{\prime }}^{*})$, i = 1, …, 105, are computed by (23), its ML tree is T$\widehat{i}$^* with ${\widehat{i}}^{*}={arg max}_{i=1,\dots ,105}{\ell }_i({\widehat{\theta }}_i^{*};{{X}}_{n^{\prime }}^{*})$.

The non-parametric bootstrap probability of tree Ti is obtained as follows. We generate B bootstrap replicates $X_{n^{\prime }}^{*b}$, b = 1, …, B. In this paper, we used B = 10⁵. For each $X_{n^{\prime }}^{*b}$, the ML tree T$\widehat{i}$^*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

$\begin{array}{ll}\text{BP}(\text{T}i,n^{\prime })=\#\{{\widehat{i}}^{*b}=i,b=1,\dots ,B\}/B.& (24)\end{array}$

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^{\prime }}^{*})$ mentioned in section 3.4 is

$\begin{array}{ll}{Y}^{*}=V_n^{-1/2}{L}_{n^{\prime }}^{*},& (25)\end{array}$

where $L_{n^{\prime }}^{*}=(1/n^{\prime }){({\ell }_1^{*},\dots ,{\ell }_{105}^{*})}^T$ with ${\ell }_i^{*}={\ell }_i({\widehat{\theta }}_i^{*};{{X}}_{n^{\prime }}^{*})$ 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 σ² = n/n′. It also follows from the above argument that $\mathrm{\text{var}}({\ell }_i^{*}-{\ell }_j^{*})\approx (n^{\prime }/n)\Vert {\xi }_i-{\xi }_j{\Vert }^2$, and thus the variance of log-likelihood difference is

$\begin{array}{ll}\text{var}\left({\ell }_i({\widehat{\theta }}_i;{{X}}_n)-{\ell }_j({\widehat{\theta }}_j;{{X}}_n)\right)\approx \Vert {\xi }_i-{\xi }_j\Vert {}^2,& (26)\end{array}$

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;{\widehat{\theta }}_i)$, $p_j(x;{\widehat{\theta }}_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).

6.2. Visualization of Probability Models

For representing the probability distribution of tree Ti, we define ${\xi }_i:={({\xi }_{1i},\dots ,{\xi }_{ni})}^T\in {ℝ}^n$ from (22) for i = 1, …, 15. The idea behind the visualization of Figure 3 is that locations of ξ_i in ℝⁿ will represent locations of $p_i(x;{\widehat{\theta }}_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)\Vert {\xi }_i-{\xi }_j{\Vert }^2$, the squared distance in ℝⁿ approximates n times Jeffreys divergence

$\begin{array}{ll}{\Vert {\xi }_i-{\xi }_j\Vert }^2\approx n\times \left(D_{\text{KL}}(p_i(x;{\widehat{\theta }}_i)\Vert p_j(x;{\widehat{\theta }}_j))+D_{\text{KL}}(p_j(x;{\widehat{\theta }}_j)\Vert p_i(x;{\widehat{\theta }}_i)\right)& (27)\end{array}$

for non-nested models (Shimodaira, 2001, section 6). When a model p₀ is nested in p_i, it becomes $\Vert {\xi }_i-{\xi }_0{\Vert }^2\approx 2n\times D_{\mathrm{\text{KL}}}(p_i(x;{\widehat{\theta }}_i)\Vert p_0(x;{\widehat{\theta }}_0))\approx 2\times ({\ell }_i({\widehat{\theta }}_i;{{X}}_n)-{\ell }_0({\widehat{\theta }}_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.

FIGURE 7

Figure 7. Three versions the visualization of probability distributions for the best 15 trees drawn using different sets of models. (A) Only the 15 bifurcating trees. (B) 15 bifurcating trees + 10 partially resolved trees + 1 star topology. This is the same plot as Figure 3B. (C) 15 bifurcating trees + 1 star topology. Note that (B,C) are superimposed, since their plots are almost indistinguishable.

For dimensionality reduction, we have to specify the origin c ∈ ℝⁿ and consider vectors a_i: = ξ_i−c. A naive choice would be the average $c=\sum _{i=1}^{15}{\xi }_i/15$. By applying PCA without centering and scaling (e.g., prcomp with option center=FALSE, scale=FALSE in R) to the matrix (a₁, …, a₁₅), 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 = η₀ for computing a_i = ξ_i − η₀ and b_i = η_i − η₀. There is hierarchy of models: η₀ is the submodel nested in all the other models, and η₁, η₂, η₃, for example, are submodels of ξ₁ (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 = β₁x₁ + ··· + β₂₅x₂₅ + ϵ of multiple regression from submodels y = β₁x₁ + ϵ, …, y = β₂₅x₂₅ + ϵ. 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,\dots ,b_{25})\in {ℝ}^{n\times 25}$ and $d={(||b_1|{|}^2,\dots ,||b_{25}|{|}^2)}^T\in {ℝ}^{25}$, then the vector for the full model is computed approximately by

$\begin{array}{ll}{a}_X=B{(B^{T}B)}^{-1}d.& (28)\end{array}$

For the visualization of the best 15 trees, we may use only b₁, …, b₁₁, 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₁, …, a₁₅ 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₁, …, a₁₅) and $d={(\Vert a_1{\Vert }^2,\dots ,\Vert a_{15}{\Vert }^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.

6.3. Asymptotic Theory of Smooth Surfaces

For expressing the shape of the region ${R}\subset {ℝ}^{m+1}$, we use a local coordinate system (u, v) ∈ ℝ^{m + 1} with u ∈ ℝ^m, v ∈ ℝ. In a neighborhood of y, the region is expressed as

$\begin{array}{ll}{R}=\left\{(u,v)\mid v\le -h(u),u\in {ℝ}^m\right\},& (29)\end{array}$

where h is a smooth function; see Shimodaira (2008) for the theory of non-smooth surfaces. The boundary surface $\partial {R}$ is expressed as v = −h(u), u ∈ ℝ^m. We can choose the coordinates so that y = (0, β₀) (i.e., u = (0, …, 0) and v = β₀), and h(0) = 0, ∂h/∂u_i|₀ = 0, i = 1, …, m. The projection now becomes the origin $\widehat{\mathit{\mu }}=(\mathbf{0},0)$, and the signed distance is β₀. The mean curvature of surface $\partial {R}$ at $\widehat{\mathit{\mu }}$ is now defined as

$\begin{array}{ll}{\beta }_1=\frac{1}{2}{\displaystyle \sum _{i=1}^m{\frac{{\partial }^{2}h\left(u\right)}{\partial u_i\partial u_i}|}_0,}& (30)\end{array}$

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 β₁ ≥ 0, whereas concave ${R}$ leads to β₁ ≤ 0. In particular, β₁ = 0 when $\partial {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 $\partial {R}$ approaches flat as n → ∞. More specifically, the magnitude of mean curvature is of order ${\beta }_1=O_p(n^{-1/2})$. The magnitude of ${\partial }^{3}h/\partial u_i\partial u_j\partial 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.

6.4. Bridging the Problem of Regions to the Z-Test

Here we explain the problem of regions in terms of the z-test by bridging the multivariate problem of section 3 to the 1-dimensional case of section 1.

Ideal p-values are uniformly distributed over p ∈ (0, 1) when the null hypothesis holds. In fact, $\mathrm{\text{AU}}({R}|Y)~U(0,1)$ for $\mathit{\mu }\in \partial {R}$ as indicated in (11). The statistic $\mathrm{\text{AU}}({R}|Y)$ may be called pivotal in the sense that the distribution does not change when $\mathit{\mu }\in \partial {R}$ moves on the surface. Here we ignore the error of $O_p(n^{-1})$, and consider only the second order asymptotic accuracy. From (10), we can write $\mathrm{\text{AU}}({R}|Y)\simeq \bar{\Phi }({\beta }_0(Y)-{\beta }_1(Y))$, where the notation such as β₀(Y) and β₁(Y) indicates the dependency on Y. Since β₁(Y) ≃ β₁(y) = β₁, we treat β₁(Y) as a constant. Now we get the normal pivotal quantity (Efron, 1985) as ${\bar{\Phi }}^{-1}(\mathrm{\text{AU}}({R}|Y))={\beta }_0(Y)-{\beta }_1~N(0,1)$ for $\mathit{\mu }\in \partial {R}$. More generally, it becomes

$\begin{array}{ll}{\beta }_0(Y)-{\beta }_1~N({\beta }_0(\mu ),1),\mu \in {ℝ}^{m+1}.& (31)\end{array}$

Let us look at the z-test in section 1, and consider substitutions:

$\begin{array}{ll}Z={\beta }_0(Y)-{\beta }_1,\theta ={\beta }_0(\mu ),c=-{\beta }_1.& (32)\end{array}$

The 1-dimensional model (1) is now equivalent to (31). The null hypothesis is also equivalent: $\theta \le 0\iff {\beta }_0(\mathit{\mu })\le 0\iff \mathit{\mu }\in {R}$. We can easily verify that AU corresponds to p(z), because $p(z)=\bar{\Phi }(z)=\bar{\Phi }({\beta }_0(y)-{\beta }_1)\simeq \mathrm{\text{AU}}({R}|y)$, which is expected from the way we obtained (31) above. Furthermore, we can derive SI from p(z, c). First verify that the selection event is equivalent: $Z>c\iff {\beta }_0(Y)-{\beta }_1>-{\beta }_1\iff {\beta }_0(Y)>0\iff Y\in {{R}}^C$. Finally, we obtain SI as $p(z,c)=p(z)/\bar{\Phi }(c)\simeq \bar{\Phi }({\beta }_0(y)-{\beta }_1)/\bar{\Phi }(-{\beta }_1)\simeq \mathrm{\text{SI}}({R}|y)$.

6.5. Model Fitting in Multiscale Bootstrap

We have used thirteen σ² 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 ${\beta }_0+{\beta }_1{\sigma }^2$ in Figure 6. Therefore we fit other models to the observed values of ${\psi }_{{\sigma }^2}$ as implemented in scaleboot package (Shimodaira, 2008). For example, poly.k model is $\sum _{i=0}^{k-1}{\beta }_i{\sigma }^{2i}$, and sing.3 model is ${\beta }_0+{\beta }_1{\sigma }^2{(1+{\beta }_2(\sigma -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 β₀ and β₁ are estimated from the tangent line to the fitted curve of ${\psi }_{{\sigma }^2}$ at σ² = 1. In Figure 6, the tangent line is drawn as red line for extrapolating ${\psi }_{{\sigma }^2}$ to σ² = −1. Shimodaira (2008) and Terada and Shimodaira (2017) considered the Taylor expansion of ${\psi }_{{\sigma }^2}$ at σ² = 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 σ² values (ten σ⁻² values: 0.5, 0.6, …, 1.4). Only the linear model ${\beta }_0+{\beta }_1{\sigma }^2$ is fitted there. The estimated β₀ and β₁ should be very close to those estimated from the tangent line described above. An advantage of using wider range of σ² in scaleboot is that the standard error of β₀ and β₁ will become smaller.

6.6. General Formula of Selective Inference

Let ${H},{S}\subset {ℝ}^{m+1}$ be regions for the null hypothesis and the selection event, respectively. We would like to test the null hypothesis $H_0:\mathit{\mu }\in {H}$ against the alternative $H_1:\mathit{\mu }\in {{H}}^C$ conditioned on the selection event $y\in {S}$. We have considered the outside mode ${H}={R},{S}={{R}}^C$ in (18) and the inside mode ${H}={{R}}^C,{S}={R}$ in (20). For a general case of ${H},{S}$, Terada and Shimodaira (2017) gave a formula of approximately unbiased p-value of selective inference as

$\begin{array}{ll}\text{SI}({H}|{S},y)=\frac{\bar{\Phi }({\beta }_0^{{H}}-{\beta }_1^{{H}})}{\bar{\Phi }({\beta }_0^{{S}}+{\beta }_0^{{H}}-{\beta }_1^{{H}})},& (33)\end{array}$

where geometric quantities β₀, β₁ are defined for the regions ${H},{S}$. We assumed that ${H}$ and ${{S}}^C$ are expressed as (29), and two surfaces $\partial {H},\partial {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 $\partial {H}=\partial {R}$ and $\partial {S}=\partial {{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}\subset {{H}}^C$, because a reasonable test would not reject H₀ unless $y\in {{H}}^C$. Note that $y\in {S}\subset {{H}}^C$ implies $0\le -{\beta }_0^{{S}}\le {\beta }_0^{{H}}$. Therefore, ${\beta }_0^{{H}}+{\beta }_0^{{S}}\ge 0$ leads to

$\begin{array}{ll}\text{SI}({H}|{S},y)\ge \text{SI}({H}|y),& (34)\end{array}$

where $\mathrm{\text{SI}}({H}|y):=\mathrm{\text{SI}}({H}|{{H}}^C,y)$ is obtained from (33) by letting ${\beta }_0^{{H}}+{\beta }_0^{{S}}=0$ for ${S}={{H}}^C$. The p-value $\mathrm{\text{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}\subset {{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}}_{\mathrm{\text{E2}}}^C$ and ${S}={{R}}_{\mathrm{\text{T1}}}$. In this case, the two surfaces $\partial {H},\partial {S}$ may not be very parallel to each other, thus violating the assumption of $\mathrm{\text{SI}}({H}|{S},y)$, so we only intend to show the potential difference between the two p-values. The geometric quantities are ${\beta }_0^{{H}}=-{\beta }_0^{\mathrm{\text{E2}}}=1.59$, ${\beta }_1^{{H}}=-{\beta }_1^{\mathrm{\text{E2}}}=-0.12$, ${\beta }_0^{{S}}={\beta }_0^{\mathrm{\text{T1}}}=-0.41$; the p-values are calculated using more decimal places than shown. SI of E2 conditioned on selecting T1 is

$\text{SI}({H}|{S},y)=\frac{\bar{\Phi }(1.59+0.12)}{\bar{\Phi }(-0.41+1.59+0.21)}=0.448,$

and it is very different from SI of E2 conditioned on selecting E2

$\text{SI}({H}|y)=\frac{\bar{\Phi }(1.59+0.12)}{\bar{\Phi }(0.12)}=0.097,$

where ${\mathrm{\text{SI}}}^{\prime }({{R}}_{\mathrm{\text{E2}}}^C|y)=1-\mathrm{\text{SI}}({{R}}_{\mathrm{\text{E2}}}^C|y)=0.903$ is shown in Table 2. As you see, $\mathrm{\text{SI}}({H}|y)$ is easier to reject H₀ than $\mathrm{\text{SI}}({H}|{S},y)$.

6.7. 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_{\mathrm{\text{all}}}=(2N-5)!/(2^{N-3}(N-3)!)$ for trees and $K_{\mathrm{\text{all}}}=2^{N-1}-(N+1)$ for edges. Next consider the number of selected regions K_select. In inside mode, regions with $y\in {{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\notin {{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 $\mathit{\mu }\notin {{R}}_i$ in inside mode and $\mathit{\mu }\in {{R}}_i$ in outside mode, and thus K_true is the same as the number of regions with $y\notin {{R}}_i$ in inside mode and $y\in {{R}}_i$ in outside mode (These numbers do not depend on the value of y by ignoring the case of $y\in \partial {{R}}_i$). Therefore, K_true = K_all − K_select for both cases.

6.8. 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 prostate-specific 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. $\widehat{M}$ is the set of selected variables, and $\widehat{s}$ represents the signs of the selected regression coefficients. We suppose that regression responses are distributed as $Y~N(\mathit{\mu },{\tau }^2{I}_n)$ where μ ∈ ℝⁿ and τ > 0. Let e_i be the ith residual. Resampling the scaled residuals σe_i (i = 1, …, n) with several values of scale σ², 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:

$\beta ={(X^{T}X)}^{-1}{X}^T\mu ,$

where X ∈ ℝ^n×p is the design matrix. We compute four types of intervals with confidence level 1 − α = 0.95 for selected variable j. $\left[L_j^{\mathrm{\text{ordinary}}},U_j^{\mathrm{\text{ordinary}}}\right]$ is the non-selective confidence interval obtained via t-distribution. $\left[L_j^{\mathrm{\text{model}}},U_j^{\mathrm{\text{model}}}\right]$ 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). By extending the method of $\left[L_j^{\mathrm{\text{model}}},U_j^{\mathrm{\text{model}}}\right]$, we also computed $\left[L_j^{\mathrm{\text{variable}}},U_j^{\mathrm{\text{variable}}}\right]$, 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

$\begin{array}{l}P({\beta }_j\in [L_j^{\text{ordinary}},U_j^{\text{ordinary}}])=1-\alpha ,\\ P({\beta }_j\in [L_j^{\text{model}},U_j^{\text{model}}]\mid \widehat{M},\widehat{s})=1-\alpha ,\\ P({\beta }_j\in [L_j^{\text{variable}},U_j^{\text{variable}}]\mid j\in \widehat{M},{\widehat{s}}_j)=1-\alpha .\end{array}$

Note that the selection event of variable j, i.e., $\{j\in \widehat{M},{\widehat{s}}_j\}$ can be represented as a union of polyhedra on ℝⁿ, and thus, according to the polyhedral lemma (Lee et al., 2016; Tibshirani et al., 2016), we can compute a valid confidence interval $\left[L_j^{\mathrm{\text{variable}}},U_j^{\mathrm{\text{variable}}}\right]$. 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 $[{\widehat{L}}_j^{\mathrm{\text{variable}}},{\widehat{U}}_j^{\mathrm{\text{variable}}}]$ by the multiscale bootstrap method of section 4 with much faster computation even for larger p.

We set λ = 10 as the penalty parameter of lasso, and the following model and signs were selected:

$\begin{array}{l}\widehat{M}=\left\{\text{lcavol},\text{lweight},\text{lbph},\text{svi},\text{pgg45}\right\},\\ \widehat{s}=(+,+,+,+,+).\end{array}$

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_j^{\mathrm{\text{model}}},U_j^{\mathrm{\text{model}}}]$ and $[L_j^{\mathrm{\text{variable}}},U_j^{\mathrm{\text{variable}}}]$, the latter is shorter, and would be preferable. This is because the selection event of the latter is less restrictive as $\{\widehat{M},\widehat{s}\}\subseteq \{j\in \widehat{M},{\widehat{s}}_j\}$; see section 6.6 for the reason why larger selection event is better. Finally, we verify that $[{\widehat{L}}_j^{\mathrm{\text{variable}}},{Û}_j^{\mathrm{\text{variable}}}]$ approximates $[L_j^{\mathrm{\text{variable}}},U_j^{\mathrm{\text{variable}}}]$ very well.

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).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors appreciate the feedback from the audience of seminar talk of HS at Department of Statistics, Stanford University. The authors are grateful to Masami Hasegawa for his insightful comments on phylogenetic analysis of mammal species.

References

Adachi, J., and Hasegawa, M. (1996). Model of amino acid substitution in proteins encoded by mitochondrial DNA. J. Mol. Evol. 42, 459–468. doi: 10.1007/BF02498640

PubMed Abstract | CrossRef Full Text | Google Scholar

Akaike, H. (1974). A new look at the statistical model identification. IEEE Trans. Autom. Cont. 19, 716–723. doi: 10.1109/TAC.1974.1100705

CrossRef Full Text | Google Scholar

Amari, S.-I., and Nagaoka, H. (2007). Methods of Information Geometry, Translations of Mathematical Monographs, Vol. 191. Providence, RI: American Mathematical Society.

Burnham, K. P., and Anderson, D. R. (2002). Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach, 2nd Edn. New York, NY: Springer-Verlag.

Cox, D. R. (1962). Further results on tests of separate families of hypotheses. J. R. Stat. Soc. Ser. B (Methodol.). 24, 406–424. doi: 10.1111/j.2517-6161.1962.tb00468.x

CrossRef Full Text | Google Scholar

Efron, B. (1979). Bootstrap methods: another look at the jackknife. Ann. Stat. 7, 1–26. doi: 10.1214/aos/1176344552

CrossRef Full Text

Efron, B. (1984). Comparing non-nested linear models. J. Am. Stat. Assoc. 79, 791–803. doi: 10.1080/01621459.1984.10477096

CrossRef Full Text | Google Scholar

Efron, B. (1985). Bootstrap confidence intervals for a class of parametric problems. Biometrika 72, 45–58. doi: 10.1093/biomet/72.1.45

CrossRef Full Text | Google Scholar

Efron, B., Halloran, E., and Holmes, S. (1996). Bootstrap confidence levels for phylogenetic trees. Proc. Natl. Acad. Sci. U.S.A. 93, 13429–13434. doi: 10.1073/pnas.93.23.13429

CrossRef Full Text | Google Scholar

Efron, B., and Tibshirani, R. (1998). The problem of regions. Ann. Sta. 26, 1687–1718. doi: 10.1214/aos/1024691353

CrossRef Full Text | Google Scholar

Felsenstein, J. (1981). Evolutionary trees from DNA sequences: a maximum likelihood approach. J. Mol. Evol. 17, 368–376. doi: 10.1007/BF01734359

PubMed Abstract | CrossRef Full Text | Google Scholar

Felsenstein, J. (1985). Confidence limits on phylogenies: an approach using the bootstrap. Evolution 39, 783–791. doi: 10.1111/j.1558-5646.1985.tb00420.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Fithian, W., Sun, D., and Taylor, J. (2014). Optimal inference after model selection. arXiv:1410.2597.

Google Scholar

Graur, D., Duret, L., and Gouy, M. (1996). Phylogenetic position of the order lagomorpha (rabbits, hares and allies). Nature 379:333. doi: 10.1038/379333a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Halanych, K. M. (1998). Lagomorphs misplaced by more characters and fewer taxa. Syst. Biol. 47, 138–146. doi: 10.1080/106351598261085

PubMed Abstract | CrossRef Full Text | Google Scholar

Halvorsen, K. (2015). ElemStatLearn: data sets, functions and examples from the book: “the elements of statistical learning, data mining, inference, and prediction” by trevor hastie, robert tibshirani and jerome friedman. R package. Available online at: https://CRAN.R-project.org/package=ElemStatLearn

Kishino, H., and Hasegawa, M. (1989). Evaluation of the maximum likelihood estimate of the evolutionary tree topologies from DNA sequence data, and the branching order in Hominoidea. J. Mol. Evol. 29, 170–179. doi: 10.1007/BF02100115

PubMed Abstract | CrossRef Full Text | Google Scholar

Kishino, H., Miyata, T., and Hasegawa, M. (1990). Maximum likelihood inference of protein phylogeny and the origin of chloroplasts. J. Mol. Evol. 30, 151–160. doi: 10.1007/BF02109483

CrossRef Full Text | Google Scholar

Konishi, S., and Kitagawa, G. (2008). Information Criteria and Statistical Modeling. New York, NY: Springer Science & Business Media. doi: 10.1007/978-0-387-71887-3

CrossRef Full Text | Google Scholar

Lee, J. D., Sun, D. L., Sun, Y., and Taylor, J. E. (2016). Exact post-selection inference, with application to the lasso. Ann. Stat. 44, 907–927. doi: 10.1214/15-AOS1371

CrossRef Full Text | Google Scholar

Linhart, H. (1988). A test whether two AIC's differ significantly. South Afr. Stat. J. 22, 153–161.

Google Scholar

Novacek, M. J. (1992). Mammalian phytogeny: shaking the tree. Nature 356, 121–125. doi: 10.1038/356121a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Posada, D., and Buckley, T. R. (2004). Model selection and model averaging in phylogenetics: advantages of Akaike information criterion and bayesian approaches over likelihood ratio tests. Syst. Biol. 53, 793–808. doi: 10.1080/10635150490522304

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosenthal, R. (1979). The file drawer problem and tolerance for null results. Psychol. Bull. 86, 638–641. doi: 10.1037/0033-2909.86.3.638

CrossRef Full Text | Google Scholar

Schennach, S. M., and Wilhelm, D. (2017). A simple parametric model selection test. J. Am. Stat. Assoc. 112, 1663–1674. doi: 10.1080/01621459.2016.1224716

CrossRef Full Text | Google Scholar

Shimodaira, H. (1997). Assessing the error probability of the model selection test. Ann. Inst. Stat. Math. 49, 395–410. doi: 10.1023/A:1003140609666

CrossRef Full Text | Google Scholar

Shimodaira, H. (1998). An application of multiple comparison techniques to model selection. Ann. Inst. Stat. Math. 50, 1–13. doi: 10.1023/A:1003483128844

CrossRef Full Text | Google Scholar

Shimodaira, H. (2000). Improving predictive inference under covariate shift by weighting the log-likelihood function. J. Stat. Plan. Inference 90, 227–244. doi: 10.1016/S0378-3758(00)00115-4

CrossRef Full Text | Google Scholar

Shimodaira, H. (2001). Multiple comparisons of log-likelihoods and combining nonnested models with applications to phylogenetic tree selection. Commun. Stat. Theory Methods 30, 1751–1772. doi: 10.1081/STA-100105696

CrossRef Full Text | Google Scholar

Shimodaira, H. (2002). An approximately unbiased test of phylogenetic tree selection. Syst. Biol. 51, 492–508. doi: 10.1080/10635150290069913

PubMed Abstract | CrossRef Full Text | Google Scholar

Shimodaira, H. (2004). Approximately unbiased tests of regions using multistep-multiscale bootstrap resampling. Ann. Stat. 32, 2616–2641. doi: 10.1214/009053604000000823

CrossRef Full Text | Google Scholar

Shimodaira, H. (2008). Testing regions with nonsmooth boundaries via multiscale bootstrap. J. Stat. Plan. Inference 138, 1227–1241. doi: 10.1016/j.jspi.2007.04.001

CrossRef Full Text | Google Scholar

Shimodaira, H. (2019). Scaleboot: Approximately Unbiased p-Values via Multiscale Bootstrap. R package version 1.0-0. Available online at: https://CRAN.R-project.org/package=scaleboot

Google Scholar

Shimodaira, H., and Hasegawa, M. (1999). Multiple comparisons of log-likelihoods with applications to phylogenetic inference. Mol. Biol. Evol. 16, 1114–1116. doi: 10.1093/oxfordjournals.molbev.a026201

CrossRef Full Text | Google Scholar

Shimodaira, H., and Hasegawa, M. (2001). CONSEL: for assessing the confidence of phylogenetic tree selection. Bioinformatics 17, 1246–1247. doi: 10.1093/bioinformatics/17.12.1246

PubMed Abstract | CrossRef Full Text | Google Scholar

Shimodaira, H., and Hasegawa, M. (2005). “Assessing the uncertainty in phylogenetic inference,” in Statistical Methods in Molecular Evolution, Statistics for Biology and Health, ed R. Nielsen (New York, NY: Springer), 463–493. doi: 10.1007/0-387-27733-1_17

CrossRef Full Text | Google Scholar

Shimodaira, H., and Maeda, H. (2018). An information criterion for model selection with missing data via complete-data divergence. Ann. Inst. Stat. Math. 70, 421–438. doi: 10.1007/s10463-016-0592-7

CrossRef Full Text | Google Scholar

Stamey, T., Kabalin, J., Johnstone, I., Freiha, F., Redwine, E., and Yang, N. (1989). Prostate specific antigen in the diagnosis and treatment of adenocarcinoma of the prostate II. Radical prostatectomy treted patients. J. Urol. 16, 1076–1083. doi: 10.1016/S0022-5347(17)41175-X

CrossRef Full Text | Google Scholar

Suzuki, R., and Shimodaira, H. (2006). Pvclust: an R package for assessing the uncertainty in hierarchical clustering. Bioinformatics 22, 1540–1542. doi: 10.1093/bioinformatics/btl117

PubMed Abstract | CrossRef Full Text | Google Scholar

Taylor, J., and Tibshirani, R. (2015). Statistical learning and selective inference. Proc. Natl. Acad. Sci. U.S.A. 112, 7629–7634. doi: 10.1073/pnas.1507583112

PubMed Abstract | CrossRef Full Text | Google Scholar

Terada, Y., and Shimodaira, H. (2017). Selective inference for the problem of regions via multiscale bootstrap. arXiv:1711.00949.

Google Scholar

Tian, X., and Taylor, J. (2018). Selective inference with a randomized response. Ann. Stat. 46, 679–710. doi: 10.1214/17-AOS1564

CrossRef Full Text | Google Scholar

Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. J. R. Stat. Soc. Ser. B (Methodol.). 58, 267–288. doi: 10.1111/j.2517-6161.1996.tb02080.x

CrossRef Full Text | Google Scholar

Tibshirani, R., Taylor, J., Lockhart, R., and Tibshirani, R. (2016). Exact post-selection inference for sequential regression procedures. J. Amer. Stat. Assoc. 111, 600–620. doi: 10.1080/01621459.2015.1108848

CrossRef Full Text | Google Scholar

Tibshirani, R., Tibshirani, R., Taylor, J., Loftus, J., and Reid, S. (2017). SelectiveInference: Tools for Post-Selection Inference. R package version 1.2.4. Available online at: https://CRAN.R-project.org/package=selectiveInference

Vuong, Q. H. (1989). Likelihood ratio tests for model selection and non-nested hypotheses. Econometrica 57, 307–333. doi: 10.2307/1912557

CrossRef Full Text | Google Scholar

Yang, Z. (1996). Among-site rate variation and its impact on phylogenetic analyses. Trends Ecol. Evol. 11, 367–372. doi: 10.1016/0169-5347(96)10041-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Z. (1997). PAML: a program package for phylogenetic analysis by maximum likelihood. Comput. Appl. Biosci. 13, 555–556. doi: 10.1093/bioinformatics/13.5.555

PubMed Abstract | CrossRef Full Text | Google Scholar