Parametric bootstrap inference for stratified models with high-dimensional nuisance specifications

Inference about a scalar parameter of interest typically relies on the asymptotic normality of common likelihood pivots, such as the signed likelihood root, the score and Wald statistics. Nevertheless, the resulting inferential procedures are known to perform poorly when the dimension of the nuisance parameter is large relative to the sample size and when the information about the parameters is limited. In many such cases, the use of asymptotic normality of analytical modifications of the signed likelihood root is known to recover inferential performance. It is proved here that parametric bootstrap of standard likelihood pivots results in as accurate inferences as analytical modifications of the signed likelihood root do in stratified models with stratum specific nuisance parameters. We focus on the challenging case where the number of strata increases as fast or faster than the stratum samples size. It is also shown that this equivalence holds regardless of whether constrained or unconstrained bootstrap is used. This is in contrast to when the number of strata is fixed or increases slower than the stratum sample size, where we show that constrained bootstrap corrects inference to a higher order than unconstrained bootstrap. Simulation experiments support the theoretical findings and demonstrate the excellent performance of bootstrap in extreme scenarios.


Introduction
Standard likelihood inference about a scalar parameter of interest is based on the asymptotic normality of likelihood pivots, such as the signed likelihood root, score and Wald statistics. This asymptotic approximation can be quite inaccurate in the presence of many nuisance parameters. An alternative, which guarantees higher accuracy, is based on the asymptotic normality of analytical modifications of the signed likelihood root, generally termed modified signed likelihood root (see for instance Severini, 2000, Chapter 7). In a two-index stratified asymptotic setting, in which both the dimension of the data and the number of nuisance parameters grow, the modified signed likelihood root has been proved to be highly accurate even in rather extreme scenarios with many nuisance parameters and very limited information (Sartori, 2003).
Parametric bootstrap methods provide an alternative assessment of tail probabilities for likelihood pivots and, in standard asymptotic settings, where the number of nuisance parameters is fixed and regularity conditions are satisfied (Severini, 2000, Section 3.4), have been shown to guarantee an equivalent level of asymptotic accuracy as analytical modifications of the signed likelihood root (see Young and Smith, 2005, Chapter 11). In particular, the two main variants of parametric bootstrap are constrained and unconstrained bootstrap (also know as conventional bootstrap). In the latter, the sampling distribution of the statistic is computed at the full maximum likelihood estimate, and in the former at the constrained maximum likelihood estimate for a given value of the parameter of interest. In standard asymptotic settings, constrained bootstrap (DiCiccio et al., 2001;Lee and Young, 2005) corrects inference about a scalar parameter in the presence of nuisance parameters to a higher order than unconstrained bootstrap. On the other hand, numerical differences are rarely detectable. Although bootstrap methods are, typically, more computationally demanding than analytical approximations to the distribution of pivots, they are available in some non-regular cases in which the modified signed likelihood root is not computable.
We investigate the properties of parametric bootstrap in models for stratified data in a twoindex asymptotic setting, where both the number q of strata and the sample size m of each stratum grow. In this setting, the usual likelihood pivots are asymptotically standard normal provided q = o(m), while the condition for the modified signed likelihood root is q = o(m 3 ) (Sartori, 2003). If q = O(m α ), then for 0 ≤ α < 1 the asymptotic normality of standard likelihood pivots still holds, with error of order O p (m (α−1)/2 ) (Sartori, 2003, formula (8)), while asymptotic normality fails in the highly stratified case with α ≥ 1. In that case, the aim of higher-order solutions is to recover first-order validity of inferential procedures.
We show here that parametric bootstrap provides valid inference when q = O(m α ), provided that α < 3. In particular, if 0 ≤ α < 1, constrained bootstrap is theoretically more accurate than unconstrained bootstrap, and both improve over standard first-order asymptotic results. On the other hand, when 1 ≤ α < 3 both variants of parametric bootstrap are equally accurate, recovering first-order accuracy with the same order of error as higher-order analytical solutions.
The theoretical results are supported by extensive simulation studies, which illustrate that parametric bootstrap is at least as accurate as use of the modified signed likelihood root, and provide evidence that constrained bootstrap can be even more accurate in some very extreme scenarios.

Background
Let l(θ; y) be the log-likelihood function for a parameter θ based on a sample y of size n, which is considered to be a realization of a random vector Y . We treat the case where the vector of parameters is partitioned as θ = (ψ, λ ) , where ψ is a scalar parameter of interest and λ is a vector of nuisance parameters, and denote byθ(y) = (ψ(y),λ(y) ) the maximum likelihood estimate of θ and byθ ψ (y) = (ψ,λ ψ (y) ) the constrained maximum likelihood estimate of θ for fixed ψ. We let U (θ; y) = ∇l(θ; y) denote the score vector, and j(θ; y) = −∇∇ l(θ; y) the observed information, with i(θ) = Eθ{j(θ; Y )} denoting the expected information. The argument θ will be dropped when no ambiguity arises, and components of vectors and blocks of matrices will be denoted by subscripts, so that for instance U ψ (θ; y) denotes the component of the score vector corresponding to ψ. Furthermore, the argument y will be dropped whenever evaluation is at the random variable Y instead of the sample y. For example, U ψ = U ψ (θ; Y ), U λ = U λ (θ; Y ), i ψψ = i ψψ (θ) and i ψλ = i ψλ (θ) are the (ψ, ψ) and (ψ, λ) blocks of i(θ), and so on.
An alternative to the asymptotic approximations to the distribution of (1), (2) and (3) is parametric bootstrap, which provides higher-order approximations for p-values, such as pr θ {R(ψ) ≤ R(ψ; y)}. There are two main variants of parametric bootstrap: i) unconstrained bootstrap where samples are simulated from the model atθ(y), and ii) constrained bootstrap where samples are simulated atθ ψ (y) (see DiCiccio et al., 2001;Lee and Young, 2005;Young and Smith, 2005, Chapter 11).
Similar results hold for S(ψ) and T (ψ) (Lee and Young, 2005;Young, 2009) with p-valuesp S 1 andp S 2 , andp T 1 andp T 2 , respectively. As Young and Smith (2005, Section 11.4) note, the theoretical advantage of constrained over unconstrained bootstrap is rarely supported by numerical evidence, because both types of bootstrap are able to equally improve over first-order results.
The advantage of bootstrap p-values in (4) and (6) over the use of analytical modifications to common statistics is that bootstrap does not require any additional, often tedious, algebraic derivation and implementation of the necessary modifications. Moreover, there are non-standard modelling settings, where R(ψ; y) is computable while R * (ψ; y) is not. One instance is when one or more components ofθ(y) are on the boundary of the parameter space. The main disadvantage of bootstrap is the additional computation that is typically required for the repeated model fits, which can be partly mitigated by parallel computing.
In some special cases, the distribution of (1), (2) and (3) depends only on ψ, so that constrained bootstrap, as well as simulating data at (ψ,λ(y) ) or even at (ψ, λ 0 ) for arbitrary nuisance vectors λ 0 , produces samples from the hypothesized model. This is the case when the model for fixed ψ is a transformation model (see Severini, 2000, Section 1.3). For instance, if y is a realization of Y = (Y 1 , . . . , Y n ) with independent and identically distributed components from a shape and scale model with generic density Hence, due to equivariance of the maximum likelihood estimator,λ and λλ 0 have the same distribution, whereλ 0 is the maximum likelihood estimator of λ based on Y 0 i 's. The same representation holds forλ ψ , so that the profile likelihood ratio has the same distribution as , which depends on ψ only. An example with a stratified gamma model is provided in the Supplementary Materials where simulation results confirm the exactness of the constrained bootstrap.

Two-index asymptotic theory for stratified models
We consider a stratified setting with q independent strata with m observations each. Therefore, the total number of observations is n = mq. The models considered here have λ = (λ 1 , . . . , λ q ) as nuisance parameter, where λ i is a stratum-specific parameter. Let y i = (y i1 , . . . , y im ) , i = 1, . . . , q, denote the vector of observations in the ith stratum and let y = (y 1 , . . . , y q ) . The vectors y 1 , . . . , y q are assumed to be realizations of independent random variables Y 1 , . . . , Y q from a parametric model with densities g 1 (y 1 ; ψ, λ 1 ), . . . , g q (y q ; ψ, λ q ), respectively. The observations within strata are also assumed to be realizations of independent random variables, so that g i (y i ; ψ, λ i ) = m j=1 g ij (y ij ; ψ, λ i ), where g ij (·) may be conditional on a covariate vector x ij . Under this specification, for fixed ψ, the likelihood has separable parameters λ 1 , . . . , λ q , so that where U i ψ is the contribution to U ψ from the ith stratum.
We work in a two-index asymptotic setting where q increases with m, as q = O(m α ), α > 0. The case α = 0 corresponds to the standard asymptotic setting. Sartori (2003, Section 4) showed that R(ψ), S(ψ) and T (ψ) are asymptotically equivalent to order o p (1) for α ≥ 0. Specifically, when 0 ≤ α < 1, the equivalence of the three quantities holds with relative error of order O p (n −1/2 ) = O p (m −(α+1)/2 ), and these are asymptotically standard normal. On the other hand, when α ≥ 1, asymptotic equivalence of R(ψ), S(ψ) and T (ψ) holds with error of order O p (m −1 ) and, more critically, the three statistics are not asymptotically standard normal, so that, for instance, Φ{R(ψ)} is not asymptotically uniform.
The derivation of the results is more straightforward for S(ψ) because the profile score is the sum of strata profile scores. However, the same results hold also for R(ψ) and T (ψ), since, as recalled above, they are both asymptotically equivalent to S(ψ). Let F θ (·) denote the distribution function of S(ψ) under θ, so that F θ (S(ψ)) is exactly uniform.
The core result of the paper is that asymptotic validity of both constrained and unconstrained bootstrap is guaranteed even in a two-index asymptotic setting, provided that α < 3, that is q = o(m 3 ). The latter condition is the same as the one required for validity of inference based on the modified signed likelihood root R * (ψ) (Sartori, 2003). In particular, we show that, when 0 < α < 1, and while, when 1 ≤ α < 3, and Hence, when 1 ≤ α < 3, the same order of error is obtained both with constrained and unconstrained bootstrap, unlike what happens with 0 ≤ α < 1. The case α = 0 corresponds to the standard asymptotic setting in which n = O(m), and (8) and (9) reduce to (7) and (5), respectively. A first intuition about why the two types of bootstrap have the same accuracy when α ≥ 1 is that the major effect of both bootstrap procedures is to remove the diverging bias term of the statistic, which overshadows any minor differences in theoretical performance that are found when 0 ≤ α < 1. A formal development of the result is given below. In the following we will concentrate on the more extreme case, i.e. α ≥ 1, while the proof of (8) and (9) for the case 0 < α < 1 is given in the Supplementary Materials. In order to prove both (10) and (11) we need some preliminary results about the distribution function F θ (x) of S(ψ) in the two-index asymptotic setting. From Sartori (2003, formula (6)), U p = U p (ψ) can be expanded as /2 ), having zero mean and variance i ψψ|λ , . Details about the orders in (12) are provided in the Appendix. When 0 ≤ α < 1 the terms in (12) are in descending order. Instead, when 1 ≤ α < 3, B becomes the leading term, followed by U ψ|λ . Finally, U ψ|λ is dominated by Re as well when α ≥ 3. In practice, when 1 ≤ α < 3, bootstrap procedures, as well as higher-order analytical solutions, are able to correct for B, so that U ψ|λ is again the leading term in expansion (12).
Let M (θ) = Eθ(S(ψ)) and Var θ (S(ψ)) be the expectation and variance of S(ψ). Asymptotic expansions detailed in the Appendix can be used to show that where For the development here, we assume that the distribution function of S(ψ) admits a valid Edgeworth expansion. Severini (2000, Sections 5.1-5.4) gives conditions and details for the extension of Edgeworth expansions for independent and identically distributed random variables to likelihood pivots, such as R(ψ), S(ψ), T (ψ). The basic requirement, in the continuous case, is that an Edgeworth expansion exists for the joint distribution of log-likelihood derivatives up to the third order, implying where the order of the remainder term is that of the third cumulant of S(ψ). Let We first focus on constrained bootstrap. From (16), The Taylor expansions in the Appendix give and where ∆ = O p (m (α−3)/2 ) is given in expression (A9) of the Appendix. Using (18) and (19), we can write As a result, if α < 3, then the following Taylor expansion of (17) holds where the error is of order O p (m −1 ) because, for α < 3, min (1, (5 − α)/2) = 1, while the error term in (17) is o p (m −1 ) whenever α > 1. In order to prove (10), note that Fθ ψ (S(ψ)) ≤ u is equivalent to S(ψ) ≤ s u , with s u the u-quantile of Fθ ψ (·), such that Fθ ψ (s u ) = u. Let s 0 u be the u-quantile of F θ (·). It is useful to express s u in terms of s 0 u . Using (20), , with ∆ given by (A9), and such that Eθ(∆) = O(m (α−3)/2 ). Moreover, we have since Var θ (∆) = O(m −2 ) and Cov θ (S(ψ), ∆) = O(m −2 ), where the order of the latter is determined by the orthogonality between U ψ|λ and the leading term of b 1 (θ) in (A7). Finally, (10) holds because The proof of (11) for unconstrained bootstrap is obtained along the same steps as above. In particular, expansion (A12) holds for Fθ(x), having the same form as (20), with ∆ replaced by ∆ 1 , which is still of order O p (m (α−3)/2 ). Details are given in the Appendix. However, while (21) is still true, (22) holds with an error of order O(m −1 ), because there is no orthogonality between U ψ|λ and the leading terms of b 2 (θ), given in (A10). Therefore, for unconstrained bootstrap we have Hence, when α ≥ 1, errors in (10) and (11) are of the same order because the O(m (α−3)/2 ) error in the mean ofS(ψ) dominates the O(m −2 ) and O(m −1 ) errors in the variance ofS(ψ) in the constrained and unconstrained cases, respectively. However, the different errors in the variance ofS(ψ) may have some effects and explain why the constrained bootstrap is sometimes numerically more accurate in extreme settings.
The arguments used in the proofs of (10) and (11) suggest that a location and scale adjustment to the statistic, as done for R(ψ) in a standard asymptotic setting by DiCiccio et al. (2001) and Stern (2006), is the key requirement to recover approximate uniformity of p-values. In this respect, a bootstrap location and scale adjustment of R(ψ), S(ψ) or T (ψ) is expected to be as effective as bootstrapping the distribution of the statistic. This conjecture is confirmed by the numerical results, both in the following section and in the Supplementary Materials.

Simulation studies
The finite-sample properties of unconstrained and constrained parametric bootstrap are assessed through extensive simulation studies, for three statistical models for stratified data. In particular, we consider a beta model, a curved exponential family model and a truncated regression model, with the results for further models reported in the Supplementary Materials. For each model, we conduct 9 simulation experiments, one for each combination of number of strata q ∈ {10, 100, 1000} and stratum sample size m ∈ {4, 8, 16}. Each simulation experiment involves 10000 simulated samples under the model at a fixed parameter vector θ 0 = (ψ 0 , λ 0 ) . For each simulated sample, 17 statistics and 6 bootstrap-based p-values are computed for testing ψ = ψ 0 . In particular, the statistics that are computed are i) R(ψ), S(ψ), T (ψ), ii) the location and location-and-scale adjusted versions of R(ψ), S(ψ), T (ψ), where the mean and variance of each statistic are estimated using unconstrained bootstrap (at θ) and constrained boostrap (atθ ψ ), and iii) R * (φ) and the signed likelihood root computed from the modified profile likelihood (see, for instance Severini, 2000, Chapter 8). The higherorder adjustment required for the latter two statistics is obtained using expected moments of likelihood quantities as in Severini (2000, Section 7.5). Finally, for each of R(ψ), S(ψ), and T (ψ), we compute the unconstrained and constrained bootstrap p-values in (4) and in (6), respectively.
In the interest of space, in what follows, we only report results for the 6 statistics based on R(ψ) shown in Table 1 The above experiments involve high-dimensional parameter spaces with as many as 1000 nuisance parameters. As a result, the assessment of the statistics requiring bootstrapping is demanding in terms of computational time and cost, even when parallel computing with a large number of cores is used. For this reason, the number of bootstrap samples is limited to 1000 in all simulation experiments.
The three blocks of rows in Table 2 give the estimated tail probabilities of the statistics of interest for the case q = 1000 and m = 8 for all three models considered. This combination of q and m was selected because it is the least extreme setting (compared to the most extreme q = 1000, m = 4) where departures from the expected behaviour in terms of the distribution of the statistics starts becoming apparent; the results for all the other combinations of q and m are provided in the Supplementary Materials. The following sections give a more detailed discussion on the figures shown in Table 2.

Beta model
As a first example, we suppose that Y ij has a beta distribution, with density function where B(·) is the beta function. The parameter of interest is ψ = log φ, whereas the stratumspecific nuisance parameters are given by The simulation experiments are carried out for ψ 0 = log(2) and the elements of λ 0 are generated from a standard normal distribution and held fixed over all the replications. The left panel of Figure 1 shows the empirical densities for the statistics in Table 1. The performance of the statistics is evaluated in terms of the closeness of their empirical density to the standard normal density. This assessment is valid also for the constrained and unconstrained bootstrap p-values, since they have been mapped into the standard normal scale by the Φ −1 (·) transformation.
The large location bias of the distribution of R(ψ) is apparent for all shown combinations of q and m, and it becomes huge for q = 1000 and m ∈ {4, 8}. All higher-order accurate statistics  Table 1 for various combinations of q and m. The N(0, 1) and Uniform(0, 1) density functions are superimposed for statistics (left) and p-values (right).
result in a marked finite-sample correction, with R * (ψ) and the unconstrained bootstrap illustrating some discrepancy from the standard normal distribution for large q/m ratios, such as q = 1000 and m ∈ {4, 8}. This is also apparent from the entries in Table 2. From the right panel of Figure 1, it is noticeable that the p-values based on R c , the location adjusted version R c l and the location-and-scale adjusted version R c ls are all close to one another. Hence, the necessary adjustment for making the distribution of R(ψ) to be close to standard normal appears to be mainly a location adjustment.

Curved exponential family
This example involves normally distributed random variables Y ij , each with mean exp(λ i ) and variance exp(ψ + λ i /2). This model was studied in Sartori et al. (1999), where it is pointed out that a marginal likelihood for ψ is not available. The simulation experiments are carried out for ψ 0 = log(1/2) and the elements of λ 0 are generated from a standard normal distribution and held fixed over all the replications.
The left panel in Figure 2 shows the empirical density functions of the statistics in Table 1, and the right panel shows the corresponding p-value distributions. Like in the previous example, the empirical, finite-sample distributions of R(ψ) are far from standard normal, while all the higher-order statistics perform considerably better. The conclusions are similar to those from the simulation experiments for the beta model, in that the required adjustment to R(ψ) seems to be a location correction. The main difference is the fact that no statistic appears to perform well for (q, m) = (1000, 4); see also the empirical tail probabilities in Table 2.  Table 1 for various combinations of q and m. The N(0, 1) and Uniform(0, 1) density functions are superimposed for statistics (left) and p-values (right).

Truncated linear regression model
The last example is taken from the econometric literature; see Greene (2004), Bartolucci et al. (2016) and the references therein. We define the response variable Y ij to be distributed as Y * ij conditionally on y * ij > 0, with y * ij = λ i + x ij ψ + ε ij , i = 1, . . . , q, j = 1, . . . , m , where the error term ε ij is standard normally distributed. For the simulation study, ψ is set to 1 and the elements of λ 0 are generated from a standard normal distribution and held fixed over all the replications. Likewise, the values x ij are generated from a standard normal distribution and held fixed over all the replications.
The left panel in Figure 3 shows the empirical density functions of the statistics in Table 1, and the right panel shows the corresponding p-value distributions. Differently from the other examples, here the distribution of the first-order statistics requires only a moderate adjustment even in the most extreme settings, and both the bootstrap-based statistics as well the R * (ψ) statistic perform rather well, providing results very close to the target distributions.

Concluding remarks
The contribution of this paper is to formally show that, in stratified settings, inference based on either unconstrained or constrained parametric bootstrap of usual likelihood pivots is effective  Table 1 for various combinations of q and m. The N(0, 1) and Uniform(0, 1) density functions are superimposed for statistics (left) and p-values (right).
in recovering their inferential performance, even in rather extreme settings, where the bias of the profile score renders vanilla first-order inference invalid. Unconstrained and constrained bootstrap for the signed likelihood ratio root, the score statistic and the Wald statistic can both recover inferential performance in stratified settings when q = O(m α ), for 0 < α < 3. As in the case for α = 0 (Lee and Young, 2005), when 0 < α < 1, constrained bootstrap is seen to have a higher degree of asymptotic accuracy than unconstrained bootstrap. On the other hand, the two bootstraps are asymptotically equivalent when 1 ≤ α < 3. The condition q = O(m α ), for 0 < α < 3, is the same as the one found in Sartori (2003) for validity of inference based on R * and on the signed likelihood root computed from the modified profile likelihood.
The results of Section 4 from the extensive simulation studies for the finite-sample assessment of the performance of constrained and unconstrained bootstrap are in par to what is expected from theory. In extreme settings, like the beta model with (q, m) = (1000, 4), constrained bootstrap appears to perform slightly better than unconstrained bootstrap. Furthermore, in all simulation experiments we carried out and as q/m diverges, the inferential performance from constrained and unconstrained bootstrap of first-order statistics seems to be deteriorating much slower than that of R * and the signed likelihood root computed from the modified profile likelihood (see also the Supplementary Materials). As a result, the evidence from the simulation studies points out that inference from parametric bootstrap is more resilient to increasing q/m than inference from well-used, analytically available higher-order statistics, with constrained bootstrap being the most accurate in extreme scenarios.
The theoretical developments in this paper do not immediately cover situations where the random variables have discrete support, because the Edgeworth expansion in (15) can only be valid for models with continuous support. The impact of discreteness on the performance of parametric bootstrap is examined in the Supplementary Materials through a binomial matched pairs model. In particular, the experimental setup of Section 4 is used for a stratified logistic regression model, where Y ij has a Bernoulli distribution with probability exp(λ i +ψx j )/{1+exp(λ i +ψx j )}, with x j = 1 for j ∈ {1, . . . , m/2} and x j = 0 for j ∈ {m/2+1, . . . , m}. The results in Figures S21-S24 and Tables S3-S11 in the Supplementary Materials indicate that the equivalence between unconstrained and constrained bootstrap of the first-order statistics in continuous models may not hold for discrete settings. In those cases, despite that unconstrained bootstrap appears to deliver a marked inferential improvement to first-order statistics, constrained bootstrap, similarly to R * , is found to perform considerably better for most combinations of q and m.
The simulation experiments in this paper have been carried out with 1000 bootstrap replications. This value is smaller than some of the recommendations of millions of replications that have appeared in the literature for standard asymptotics settings (Young, 2009;DiCiccio et al., 2017). For stratified settings with α > 1 the bootstrap adjustments have the role of recovering asymptotic uniformity of p-values, rather that providing a small-sample refinement of p-values that are asymptotically uniform. As a result, use of a huge number of bootstrap replications is less essential, and the few experiments we carried out with more than 1000 bootstrap replications are in support of that statement. More comprehensive simulation studies to support that statement are unfortunately not feasible with current computing capabilities.

Supplementary Materials
The Supplementary Materials provide the outputs from the simulation experiments described in Section 4, for all models and all combinations of statistics, q and m. Outputs are also provided for other models, given by a gamma model, a Behrens-Fisher model and by the logistic regression model described in Section 5. The outputs include null distributions of the various statistics and distributions of p-values, through extended versions of the Figures 1-3, and empirical tail probabilities, through extended versions of Table 2. random variables X 1 , . . . , X q . Then where U i ψ is the contribution to U ψ from the ith stratum, and U λ = (U λ 1 , . . . , U λq ) . Here and in the following, when the argument is omitted, evaluation at θ is understood.
The terms on the right-hand side of (12) are seen to be of order Derivation of (13) and (14) As a first step, consider the expansion where when α > 1 both terms in C are of the same order, which is again determined using (A1). Hence, with C/i ψψ|λ = O p (m −1 ).

Derivation of (11)
First, from (16), we have In order to obtain expansions of M (θ) and v(θ) around θ, we use the fact that, when α > 1, (Sartori, 2003). This implies that an expansion for Fθ(x) of the form (20) holds with a different ∆ term, which is still of order O p (m (α−3)/2 ).
In order to obtain an expansion for M (θ) we follow the same steps as in (A6)-(A9), giving (18). In particular, we have , so that the first three summands on the right hand side of the last formula are of order O p (m α−1 ), while the remaining two are of order O p (m α−2 ). This leads to where the term ∆ 1 is of order O p (m (α−3)/2 ), as its expected value, because the leading terms in (A10) are of the same order as b 1 (θ) in (A6). Using (A11) and an expansion similar to (19) we obtain so that the same error as in (20) holds also for unconstrained bootstrap, i.e.
The steps leading from (A12) to (11) are the same as those from (20) to (10).
References Barndorff-Nielsen, O. E. (1986). Inference on full or partial parameters based on the standardized signed log likelihood ratio. Supplementary Material for S1 Proof of (3.8) and (3.9) In the same setting as in Section 3 of the main text, we consider q = O(m α ), 0 ≤ α < 1. The case α = 0 corresponds to the standard asymptotic setting, where n = O(m).
The core results we prove here are Hence, as with α = 0 (Lee and Young, 2005), the order of error is different with constrained and unconstrained bootstrap.
Higher order cumulants of S(ψ), r = 3, 4, . . ., have the form as in standard asymptotics. As in Section 3 of the main text, the developments here rely on the assumption that the distribution function of S(ψ) admits a valid Edgeworth expansion. In particular, in the continuous case, where the order of the remainder term is that of the third cumulant of S(ψ). Using (S1.4) and (S1.6) we have We first focus on constrained bootstrap. From (S1.11), (S1.12) We, then, show the analogous of expression (3.18) of the main text in the case 0 ≤ α < 1, i.e.
The remaining steps are analogous to those from (S1.20) with ∆ replaced by ∆ 1 .

S2 Full set of simulation results
We report the full set of results of the simulation studies described in the main text. In addition to the models described there, which include also the logistic regression model described in the final section, two further models are considered, given by a gamma model and a Behrens-Fisher model, respectively. Some details about them are given as follows.

Gamma model
We take a stratified gamma model with common shape parameter, for which Y ij has density function where Γ(·) is the gamma function. The parameter of interest is ψ = log(α) and the stratum-specific nuisance parameters are given by λ i = log(β i ). The simulation data sets where generated for ψ 0 = log(2) and the elements of λ 0 are fixed to the logarithm of random draws from an exponential distribution with rate 1/2 and held fixed over all the replications.
The results for this model, reported in the following, are similar to those for the beta model.

Behrens-Fisher model
Suppose that Y ij are normally distributed, with common mean ψ across all strata, and stratum-specific variances exp(λ i ). This example was studied in Young (2009), where it is illustrated that R * (ψ) is performing worse than constrained bootstrapping of R(ψ) in a simulation study with q = 20 strata and the stratum sample size varying between 3 and 10. The simulation experiments are carried out for ψ 0 = 0 and λ 0 generated from a uniform distribution in (0, 1) and held fixed over all the replications.
The findings in Young (2009) are confirmed also in the more extensive simulation studies here, whose results are reported in the following. For the Behrens-Fisher model, the location-adjustment to R(ψ) does not have much of an effect to the distribution of the statistic. On the other hand, the location-and-scale adjusted statistic recovers inferential performance. So, it seems that only a scale adjustment is sufficient to recover first-order inferential performance. Indeed, for this model the profile score bias is exactly equal to zero, so that the distribution of the first-order statistic R(ψ) is essentially centred around zero. Also, both constrained and unconstrained bootstrap perform extremely well across all the simulation scenarios, and, along with R c ls , they outperform R * (ψ) in the most extreme settings.
A remarkable observation in this model is the excellent accuracy of the score statistic S(ψ) for all combination of q and m considered. This can be explained by the fact that the leading term U ψ|λ of U p in (S1.3) is exactly normally distributed. Moreover, M (θ) and v(θ), in (S1.4) and (S1.6) respectively, are equal to zero.

Results
For each model, simulation results for the 23 statistics are visualized. For plotting purposes, the statistics have been split in two groups. Table 1 lists the first group, and Table 2 the second one.

Statistic
Plotting Symbol Description R m Signed likelihood root computed from the modified profile likelihood Table 2: Statistics computed in the simulation experiments (second group for plotting). The meanμ R and the standard deviationσ R of R(ψ) are estimated through unconstrained bootstrap, by simulating from the model at θ =θ, and likewise for S(ψ) and T (ψ). The meanμ R and the standard deviationσ R of R(ψ) are estimated through constrained bootstrap, by simulating from the model at θ =θ ψ , and likewise for S(ψ) and T (ψ).

Statistic Plotting Symbol
Description  Table 1 for various combinations of q and m. The N(0,1) density function is superimposed.  Table 1 for various combinations of q and m. The Uniform(0,1) density function is superimposed.  Table 2 for various combinations of q and m. The N(0,1) density function is superimposed.  Table 2 for various combinations of q and m. The Uniform(0,1) density function is superimposed.  Table 1 for various combinations of q and m. The N(0,1) density function is superimposed.  Table 1 for various combinations of q and m. The Uniform(0,1) density function is superimposed.  Table 2 for various combinations of q and m. The N(0,1) density function is superimposed.  Table 2 for various combinations of q and m. The Uniform(0,1) density function is superimposed.  Table 1 for various combinations of q and m. The N(0,1) density function is superimposed.  Table 1 for various combinations of q and m. The Uniform(0,1) density function is superimposed.  Table 2 for various combinations of q and m. The N(0,1) density function is superimposed.  Table 2 for various combinations of q and m. The Uniform(0,1) density function is superimposed.  Table 1 for various combinations of q and m. The N(0,1) density function is superimposed.  Table 1 for various combinations of q and m. The Uniform(0,1) density function is superimposed.  Table 2 for various combinations of q and m. The N(0,1) density function is superimposed.  Table 2 for various combinations of q and m. The Uniform(0,1) density function is superimposed.  Table 1 for various combinations of q and m. The N(0,1) density function is superimposed.  Table 1 for various combinations of q and m. The Uniform(0,1) density function is superimposed.  Figure 19: Behrens-Fisher model. Estimated null distribution of statistics for the statistics in Table 2 for various combinations of q and m. The N(0,1) density function is superimposed.  Table 2 for various combinations of q and m. The Uniform(0,1) density function is superimposed.  Table 1 for various combinations of q and m. The N(0,1) density function is superimposed.  Table 1 for various combinations of q and m. The Uniform(0,1) density function is superimposed.  Table 2 for various combinations of q and m. The N(0,1) density function is superimposed.  Table 2 for various combinations of q and m. The Uniform(0,1) density function is superimposed.

S3 Empirical tail probabilities
Tables 3-11 report empirical tail probabilities for the statistics in Table 1 of the main text,  for all models and all simulation settings. The rows of each table are subdivided into 6 different blocks, each one corresponding to a different model.