Penalized Fieller's Confidence Interval for the Ratio of Bivariate Normal Means

Abstract

Constructing a confidence interval for the ratio of bivariate normal means is a classical problem in statistics. Several methods have been proposed in the literature. The Fieller method is known as an exact method, but can produce an unbounded confidence interval if the denominator of the ratio is not significantly deviated from 0; while the delta and some numeric methods are all bounded, they are only first-order correct. Motivated by a real-world problem, we propose the penalized Fieller method, which employs the same principle as the Fieller method, but adopts a penalized likelihood approach to estimate the denominator. The proposed method has a simple closed form, and can always produce a bounded confidence interval by selecting a suitable penalty parameter. Moreover, the new method is shown to be second-order correct under the bivariate normality assumption, that is, its coverage probability will converge to the nominal level faster than other bounded methods. Simulation results show that our proposed method generally outperforms the existing methods in terms of controlling the coverage probability and the confidence width and is particularly useful when the denominator does not have adequate power to reject being 0. Finally, we apply the proposed approach to the interval estimation of the median response dose in pharmacology studies to show its practical usefulness.

1 Introduction

The ratio estimate, defined as the ratio of means of two random variables, is often encountered in biomedical studies (Tin, 1965). Many problems reduce to forming a confidence interval (CI) for the ratio of two (asymptotically) normal random variables, such as noninferiority and bioequivalence assessment, and dose-response analysis (Vuorinen and Tuominen, 1994; Faraggi et al., 2003; Pigeot et al., 2003). Hence, several statistical methods have been proposed to address this problem (Fieller, 1954; Hayya et al., 1975; Choquet et al., 1999; Wang et al., 2015). Among them, the Fieller method and the delta method are two commonly used analytic approaches (Herson, 1975; Hirschberg and Lye, 2010; Sherman et al., 2011). The Fieller method, which gives an exact CI for achieving the required coverage probability (CP), is based on the inverting of the t statistic. This type of CI is asymmetric around the ratio estimate, which is a good property, as it reflects the skewness of the small-sample distribution of the ratio. However, when the denominator of the ratio is not significantly deviated from 0, the Fieller CI will be unbounded, being either the entire real line or the union of two disconnected infinite intervals (Fieller, 1954). Moreover, the Fieller method always has a positive unbounded probability (UP) even in large samples, and its UP increases when the significance level becomes smaller. That is, when the confidence level is set to be sufficiently close to 1, the Fieller CI can always be unbounded. On the other hand, the delta method directly adopts the first-order Taylor expansion by treating the ratio as a nonlinear function of the numerator and the denominator. By assuming asymptotic normality in large samples, this method produces a symmetric and bounded CI in contrast to the Fieller method. However, the delta method often has an inaccurate CP and unbalanced tail errors even in a moderate sample size. Furthermore, previous studies have shown that the delta method has even poorer performance if the true value of the ratio has opposite sign to the correlation coefficient between the numerator and denominator (Hirschberg and Lye, 2010). Because the exact CI for a ratio estimate is naturally asymmetric, this method has been regarded having the “most serious error” in CI construction due to its enforced symmetry (Efron and Tibshirani, 1994). Although a method based on second-order Taylor expansion (Hayya et al., 1975) has also been proposed, it is rarely used because this method has the same disadvantages as the delta method, and its calculation is more involved. Because of the main drawbacks of the preceding methods, some numerical procedures have been proposed, such as the direct integral method, the parametric bootstrap (PB) method, and the nonparametric bootstrap method (Choquet et al., 1999; Wang et al., 2015). Although the CIs obtained from those iterative procedures are bounded, and generally have a better CP than the delta method, they usually result in a much wider confidence width (CW) and are more computationally intensive (Wang et al., 2015). It should be noted that all the preceding methods, other than the Fieller CI, are first-order correct (Wang et al., 2015). This fact indicates the relatively limited improvement of the numeric methods over the delta method, because they are of the same order correctness.

In this article, we propose a novel analytic approach for constructing the CI for the ratio estimate. Our method, the penalized Fieller (PF) method, uses the same principle as the Fieller method, but adopts a penalized likelihood approach to estimate the denominator. The PF method always produces a bounded CI with an appropriately chosen penalty parameter. Moreover, the proposed method is second-order correct under the bivariate normality assumption. That is, its coverage probability will converge to the nominal level faster than other bounded methods. Simulation results show that the PF method generally outperforms the existing methods in terms of controlling the CP and the CW and is particularly useful when the denominator does not have adequate power to reject being 0. Finally, we apply the new approach to the interval estimation of the median response dose, which is commonly used in pharmacology, to show its practical usefulness.

2 Method

2.1 Definitions

Suppose that a pair of random variables jointly follow a bivariate normal distribution with marginal distributions, respectively, being and ⁠, and their covariance being σ₁₂. Denote to be the correlation coefficient, and let ⁠, ⁠, ⋅⋅⋅, be an independent and identically distributed sample from the above distribution. Then, the sample means, variances, covariance, and correlation coefficient can be estimated by ⁠, ⁠, ⁠, ⁠, ⁠, and ⁠, respectively. Further, the variance and covariance of and ⁠, respectively, are ⁠, ⁠, and with their estimates, respectively, being ⁠, ⁠, and ⁠. We assume that μ₂ is a nonzero constant so that can be well defined. The commonly used point estimate for r is the ratio estimate ⁠.

2.2 Overview of the Fieller Method

Because the PF method employs the same principle as the Fieller method, we first briefly introduce the latter. The Fieller CI is calculated by inverting the t statistic for testing the null hypothesis ⁠. It is shown that follows a t distribution with degrees of freedom under the null hypothesis. Hence, the confidence region of the Fieller CI for r with the confidence level is identical to the set of r that leads to the acceptance of the null hypothesis at the significance level α, and is determined by the inequality ⁠, where denotes the quantile of the t distribution with degrees of freedom. Rearranging this inequality results in the following quadratic inequality with respect to r: ⁠, where ⁠, ⁠, and ⁠. The solution of this inequality depends on the sign of A and ⁠. Through simple calculation, Δ can be further expressed as ⁠. Hence, also implies ⁠. Let and (⁠⁠) be the two roots of the equation if ⁠. In the situation of ⁠, the Fieller CI is ⁠. Notably, indicates that the denominator is significantly different from 0. However, if ⁠, the Fieller CI will be unbounded. For instance, if ⁠, the Fieller CI will be ⁠; otherwise, it will be ⁠. Meanwhile, it is clear to see that A will always be less than 0 (ie, unbounded) for a sufficiently small α. It should be emphasized that the null hypothesis is equivalent to only when ⁠. However, the Fieller method does not take this information into consideration. As such, the Fieller CI has the potential to overestimate the CW.

2.3 Penalized Estimate for the Denominator

The Fieller method can be unbounded when the denominator is not significantly different from 0. This scenario will occur if the denominator has a large variance. A natural idea to address this problem is to penalize the parameter μ₂, so that we can reduce the variance of the denominator by introducing some bias. Specifically, by utilizing the information that μ₂ cannot be 0, we define a penalized log-likelihood function of as

where is the penalty parameter. It is shown in Appendix  A that the preceding penalized log-likelihood function attains its maximum value when ⁠. Plugging into the equation, we obtain the penalized estimate for μ₂ as ⁠. Without loss of generality, we assume ⁠. It can be seen that reduces to be when ⁠, and for ⁠, is a biased estimator for μ₂, and the bias is ⁠. Through the penalty, always shrinks away from 0. Further, can be rewritten as ⁠, where is the sample coefficient of variation of Y. As such, tends to penalize the large coefficient of variation of the denominator. Once is obtained, a penalized ratio estimate for r can be proposed as ⁠.

Note that can be viewed as a bivariate function of and ⁠, denoted by ⁠. By making a Taylor expansion for around μ₂ and v₂, it is not difficult to show that and ⁠, where and ⁠. It is interesting that the error terms are only instead of the regular ⁠. Note that the magnitude of measures the extent of penalty with its range being (0.5,1). Specifically, represents the case (ie, no penalty), and corresponds to the scenario ⁠. As such, we asymptotically have ⁠, which means that the variance of the denominator is reduced through the penalty. However, for the fixed penalty parameter λ, we have ⁠. Notably, can also be written as ⁠, where is the coefficient of variation of Y. Hence, the value of only depends on λ and ⁠. Since and v₂ are unknown, we estimate as ⁠, where ⁠. Similarly, we have and ⁠.

2.4 Adjustment of the Numerator

Recall that the Fieller method is based on the inverting of the t statistic where the null hypothesis is ⁠. To construct the test statistic, an unbiased estimator is used to estimate for the Fieller's method. However, if we simply replace by ⁠, then the latter is a biased estimator for any ⁠. Further, it is shown that ⁠. Let be any of the root-n consistent estimates for r, and define ⁠. If we adjust by ⁠, it can be seen that ⁠. Hence, is expected to have a smaller bias than ⁠. Note that is an exception, since is unbiased in this situation. If ⁠, the bias of is ⁠, while the bias of is generally ⁠, which would be asymptotically smaller.

On the other hand, we need to select an appropriate point estimate so as to make the bias of as small as possible. There are many kinds of point estimates for r that we can choose for such an adjustment (Tin, 1965). For simplicity, we only consider the commonly used ratio estimate and the proposed penalized ratio estimate as the plug-in estimates in this article. It is shown that and ⁠. By the facts and ⁠, we know that the difference between and should be small in large samples, but the latter provides a more conservative adjustment. In order to comprehensively compare the bias-corrected estimates and ⁠, Web Tables 1 and 2 provide the estimated biases of ⁠, ⁠, and under various simulation configurations. From those tables, we can see that usually has the smallest bias when ⁠, and generally has less bias than when λ is set at the recommended values (see Section 2.5), especially in the case of small n and large ⁠. Based on these results, we recommend adjusting by (denoted by ⁠). Moreover, has a dramatically simple expression, which also makes it easy to further evaluate its variance. Since and ⁠, we, respectively, estimate and as and ⁠. Note that is only a term of ⁠, but we still keep it in so that both and can always be guaranteed, where is the estimated correlation coefficient for and ⁠.

2.5 Selecting a Suitable Penalty Parameter

Before computing the PF CI, we need to select an appropriate penalty parameter λ. When ⁠, the PF method reduces to the Fieller method. However, if ⁠, it is shown in Appendix  B that the CW of the PF CI will tend to be 0, and the rate of convergence is of order ⁠. This fact indicates that, when ⁠, the CP of the PF method converges to 0 as well. Therefore, when λ is too small, the PF CI will still have a positive UP. However, if we select a very large λ, then the PF CI may risk underestimating the CP. Selecting λ is a trade-off between the UP and the CP, we have the following theorem.

Theorem 1 The PF CI has a lower UP than the Fieller CI, and specially, the former is always bounded (ie, UP=0) if and only if the penalty parameter ⁠.

The proof of Theorem 1 is provided in Appendix C. From the theorem, is the minimum value that guarantees the bounded property of the PF CI. As such, it is reasonable to believe that selecting is a good choice for the trade-off between the UP and the CP. To verify this, Web Tables 3-6 provide the empirical results of the performance of the PF CI against various λ. As shown in these tables, the PF CI tends to have a narrower CW with a larger λ, but may have more chance to underestimate the CP as λ increases. For instance, despite being slightly conservative, the PF CI with can control the CP well even when the UP of the Fieller CI is as large as 40%. We also find that the PF CI with has nearly the same performance in terms of controlling the CP as the PF CI with in our simulation study. However, when the UP of Fieller CI is larger than or equal to 10%, we have observed that the CP of the former is uniformly smaller than that of the latter. This observation indicates that the latter may maintain the valid CP in more situations. Further, when λ increases to ⁠, the PF CI may have an inflated CP even when the UP of the Fieller CI is 20% as shown in Web Tables 5 and 6, despite the fact that its CW is much narrower in some scenarios. In the case of ⁠, selecting λ becomes a trade-off between the CP and the CW. Because CP might be the most crucial criterion to evaluate a CI, we would rather select ⁠, which is shown to be the most conservative choice. Note that with the pre-set confidence level. Hence, our previous conclusions are not affected by selecting ⁠. Further, does not depend on the effect size, and can be pre-determined as long as we know the sample size n. In some circumstances, even when n is unknown, we can select as an approximation, where is the quantile of the standard normal distribution. For example, when and ⁠, λ should be 1.10. However, if we do not know the sample size, we can simply use instead.

2.6 Penalized Fieller's Confidence Interval

Once λ is determined, similar to the Fieller method, the confidence region of the PF method can be constructed by solving the following quadratic inequality: ⁠, where ⁠, ⁠, and ⁠. Note that both a and will be greater than 0 with ⁠. As such, the PF CI is directly given by ⁠.

Like other bounded methods, the PF method provides an asymptotic CI. It is known that the delta and existing numeric methods are all first-order correct. However, for the PF method, we have the following theorem.

Theorem 2 Under the bivariate normality assumption, the proposed PF method is second-order correct, that is, the differences of confidence limits between the PF method and the exact method are ⁠.

The proof of Theorem 2 is in Appendix  D. It should be emphasized that second-order correct is stronger than second-order accurate, because the former implies the latter (DiCiccio and Efron, 1996). That is, the actual CP of the PF method should be the nominal level ⁠. To further elaborate on the second-order correctness of the PF method, Web Figures 1-3 plot the estimated CP of the PF, Fieller, delta, and PB methods under the setting of (very large variation for the denominator). It is clearly shown from these figures that the CP of the PF method generally converges to the nominal level faster than those of the delta and PB methods except when and for the PB method, where the actual CP of the latter are only the nominal CP ⁠. Note that the Fieller method generally keeps an accurate CP regardless of the values of n, ⁠, and ⁠, because it is an exact method in the bivariate normality case. However, this method can have a large UP in our settings. For example, when ⁠, the UP for the Fieller method are 30.60% and 56.09% for and 0.99, respectively, which are large proportions. Hence, the PF method has an excellent performance in terms of controlling the CP and the UP.

It should be noted that the second-order correctness of the PF method depends on the assumption of bivariate normality. Under this condition, the sample means of the numerator and denominator also exactly jointly follow the bivariate normal distribution. However, if the bivariate normality assumption is violated, the sample means are only asymptotically normal. In this situation, all the methods, including the PF and Fieller methods, are first-order correct. However, it is well known that the t test is robust against nonnormality, especially when the underlying distribution is symmetric. Because the PF method is based on the inverting of the t statistic, it is reasonable to believe that the PF CI can remain a valid CP under a wide range of scenarios. We evaluate how sensitive the PF CI is to the violation of the bivariate normality assumption in Section 3. In real applications, if the numerator and the denominator are known to be independent, we simply set for computing the confidence limits of the PF method. Further, when the numerator and the denominator have different degrees of freedom, the Welch-Satterthwaite equation can be used to approximate the combined degree of freedom associated with the t statistic.

3 Simulation Study

3.1 Simulation Settings

To evaluate the proposed PF method, an extensive simulation study has been conducted for comparison with the Fieller, delta, and PB methods. For the PB method, we use a similar procedure to that in Wang et al. (2015) for calculating its confidence limits. The details of the delta and PB methods can be found in Web Appendix  A. In our simulation study, we fix the penalty parameter to be ⁠. The data are simulated from the bivariate normal distribution, the bivariate Laplace distribution, and the bivariate skew normal distribution. Notably, both the normal distribution and the Laplace distribution are symmetric, while the skew normal distribution is asymmetric. For the latter, we fix the shape parameter to be 10, so that the skewness of such distribution is about 0.96. From simulation study (data not shown), we found that the value of r has little impact on the simulation results except when ⁠. Hence, r is fixed to be 1 and 0. We set the mean of the denominator μ₂ to be 1; otherwise, we can simultaneously divide the denominator and numerator by μ₂ for any ⁠, so that the transformed denominator always has a mean of 1. For the numerator, its coefficient of variation is set to be 0.4 for ⁠, and when ⁠, we fix its variance to be 1. The sample size n is selected to be 20 and 50, and the confidence level is set to be 0.95 and 0.99. Furthermore, we have a similar finding to that of Wang et al. (2015), that is, the power of the t test (denoted by ⁠) for testing the denominator being 0 plays a prominent role in the performance of all the methods. We set to be ⁠, ⁠, 0.99, 0.90, 0.80, 0.60, 0.40, and 0.20, where such power ranges from perfect to low. Note that denotes the probability that the denominator is not statistically different from 0, which is also the theoretical UP for the Fieller method. Once n, ⁠, and are fixed, we can calculate by solving the equation ⁠, where is the distribution function of the noncentral t distribution with the noncentral parameter being ⁠. Web Table 7 gives the corresponding values of for various combinations of n, ⁠, and ⁠. It can be seen from Web Table 7 that the values of for are uniformly larger than those for when is fixed. For the choice of the correlation coefficient ρ, we emphasize the scenario in which ρ has the opposite sign to r. Therefore, ρ is assigned to be 0, −0.4, and −0.8. We also conduct a simulation study (data not shown) considering ρ being 0.4 and 0.8, but we find that the performances of all the methods are similar to that of ρ being 0. Note that when generating the data from the bivariate skew normal distribution, we first need to perform a parameterization, and the details are given in Web Appendix  B (Azzalini, 2005). Finally, the number of simulations k is set at 10 000.

We compare the performance of four types of CI based on CP, ML/(ML+MR), UP, and CW, where ML and MR, respectively, are the left and right tail errors missing the true value of r. CP is estimated by the proportion that the CI contains the true value of r among k replicates. ML and MR are calculated by and ⁠, respectively, where denotes the counting measure, and and are the confidence limits of the estimated CI. Note that means that the CI is bounded. As such, we only consider the bounded CIs when estimating the ML and the MR, as it is impossible to distinguish between the left side and the right side if the CI is the union of two disconnected infinite intervals. Further, UP is computed as ⁠. Note that the UPs of the PF, delta, and PB methods always stay at 0. On the other hand, as the Fieller CI may be unbounded (ie, has infinite length), we use the median length of CIs among k replicates to estimate the CW. It is believed that a good CI should control the CP well, have a narrow CW, and have balanced ML and MR. We deem that underestimating the CP is undesirable because of the inflated size, while slightly overestimating the CP can be acceptable. Further, if a balance between ML and MR is achieved, then ML/(ML+MR) should be close to 0.50.

3.2 Simulation Results

Tables 1 and 2, respectively, display the estimated UP, CP, CW, and ML/(ML+MR) of the two-sided 95% and 99% CIs for the PF, Fieller, delta, and PB methods when and ⁠. From these tables, it can be seen that the estimated UP of the Fieller CI is generally around the theoretical UP in all the scenarios. The Fieller method has a very accurate CP in scenarios of symmetric distributions regardless of the values of n, ⁠, and ⁠. However, if the underlying distribution is skew normal, we have observed that the Fieller method can slightly underestimate the CP when ⁠. Unlike the Fieller method, the performances of the PF and delta methods seem to be less affected by the underlying distribution. The PF CI can control the CP well as long as ⁠, regardless of the sample size, the confidence level, and the underlying distribution. However, for the delta method, we find that its CP may be underestimated when ⁠. This fact indicates that the delta method can only control the CP when the denominator has the perfect power to reject being 0 (eg, ⁠). We also observe that the CP of the PB method depends on the confidence level. When ⁠, the PB methods can control the CP, except for being a little conservative when the underlying distribution is skew normal. However, if the confidence level increases to 0.99, the PB method may underestimate the CP when ⁠, but its CP is still closer to the nominal level than that of the delta method. On the other hand, it can be seen that the delta method has the narrowest CW, while the CW of the PB method generally is the widest among all the bounded methods. The PF method typically has a width in between the delta method and PB method, and its CW is narrower than that of the Fieller CI when ⁠. In the case where the UP of the Fieller CI is larger than 50%, its median length becomes infinite as shown in Tables 1 and 2. From the values of ML/(ML+MR), we find that the delta and PB methods demonstrate the unbalanced tail errors on the left and right, because almost all the values of ML/(ML+MR) are far away from 0.50, while the PF and Fieller methods have much more balanced tail errors. This is not surprising because both the delta and PB methods are based on the normality assumption of the ratio estimate, which may not always be appropriate as it generally has a skewed distribution even in moderate samples. When the UP of Fieller's CI ranges from small to large, the tail errors of all the methods tend to be more skewed.

Notably, when ⁠, all the methods have nearly the same CP and the same CW in all the scenarios; but if ⁠, all the methods uniformly show poor performance. The interesting scenarios lie on the parameter space that is between 0.60 and 0.99. For simplicity, we restrict our following discussion to this parameter space. It can be seen that the PF method can overestimate the CP in some cases, but its CW still grows much slower than those of the Fieller and PB methods as the decreases. For the PB method, in the situation of ⁠, where it can control the CP, its CW is even wider than the Fieller CI. Hence, the PF method is superior to the Fieller and PB methods because its CI is narrower while still remaining a valid CP. Note that the delta method has a narrower CW than the PF method, but its CP is underestimated at the nominal level. When looking at 95% CI, we find that when decreases 0.10, the under-coverage of the delta CI will increase about 0.005-0.01 (10-20% inflation in type 1 error), and the CW of the PF CI is nearly 15% wider than the delta CI. When reduces from 0.99 to 0.60, the under-coverage of the delta method can range from 0.01 to 0.05 (20-100% inflation in type 1 error), despite the fact that its CW is about 15-60% narrower than that of the PF method. For the 99% CI, the delta method has an even worse CP (100-300% inflation in type 1 error), but the PF method still maintains the correct level. Recall the fact that the delta method typically has unbalanced tail errors. Hence, this method seems to narrow its width by sacrificing the accuracy of CP, and, to some extent, the accuracy of the left and right tail errors. Moreover, from Web Tables 3-6, we find that the PF method with penalty parameter demonstrates a similar CW to the delta method, but the former has a much better CP and more balanced tail errors than the latter.

3.3 Additional Simulation Results and Conclusions

All the other results of UP, CP, CW, and ML/(ML+MR) with and −0.8, or are provided in Web Tables 8-17. From Web Tables 8-11, it can be seen that ρ has little impact on the performance of the PF and Fieller methods. But when ρ changes from 0 to −0.4, and −0.8, we find that the delta and PB methods tend to have more skewed tail errors, although their CPs seem to be less affected. Web Tables 12-17 display the results for ⁠. From these tables, we find that the PF method performs even better in terms of controlling the CP when in the symmetric distribution scenarios, because it can still have an accurate CP even when ⁠. However, both the PF and Fieller methods can slightly underestimate the CP in the case of skew normal distribution. It can also be seen that the delta and PB methods tend to seriously overestimate the CP. However, the CW of the delta method is still the narrowest. The simulation results suggest that the delta method seems to be preferred only when ⁠, which, however, is just a special case.

In summary, in nearly all the simulation settings that we considered, the Fieller method generally keeps a good CP, but its CW can be overlong and even infinite; the PF method can generally control the CP well as long as the Fieller CI has a finite median CW (ie, %), while providing a competitive confidence length; although the delta method has the narrowest CW among all the methods, this may be due to its inaccurate CP and unbalanced tail errors; and the PB method generally has the worst performance among all the methods, because its CW is nearly as long as that of the Fieller method, but it can still underestimate the CP in some scenarios like the delta method. Further, the PB method is a resampling-based approach, and is more computationally intensive than other methods. Based on all the results, the PF method generally outperforms the existing methods in terms of controlling the CP and the CW and is particularly useful when the denominator does not have adequate power to reject being 0.

4 Application to the Interval Estimation of the Median Response Dose

4.1 Background

The median response dose is frequently used in pharmacology studies, especially in toxicity experiments, such as the median effective dose (⁠⁠, the dose of a drug that produces a desired effect in 50% of a tested population) or the median lethal dose (⁠⁠, the dose required to kill half the tested population after a specified duration). In such experiments, a small group of animals is exposed to a measurable stimulus at each of a small number of given doses, usually evenly distributed on a logarithmic scale. After that, the outcome of interest, usually a binary trait, such as death or affected, is recorded for each animal. For simplicity, we assume that the experiment comprises q dose levels with the base-10 log-doses l₁, l₂, …, ⁠, and for each dose level, we respectively have the replications m₁, m₂, …, ⁠. The number of the response animals at the ith level is denoted by ⁠. Consequently, s₁, s₂, …, are mutually independent binomial random variables with ⁠. We further assume that the probability is related to the log-dose through the logistic model

where represents the median response dose. For fixed γ, the slope β determines the shape of the response curve. When β ranges from 0 to ⁠, the response curve will become steeper and steeper. In the extreme cases of 0 and ⁠, will uniformly tend to be 0.5 and 1 for all the dose levels, respectively. In both situations, the data contain little information about γ. By fitting the standard logistic model, it is easy to obtain ⁠, ⁠, and their covariance matrix. The point estimate of γ is given by ⁠. It can be seen that both the numerator and denominator of come from the logistic regression coefficients. Because both and asymptotically follow normal distributions, we will show how the proposed PF method can be applied to the interval estimation of γ to three real-world datasets—the Hewlett data, the Puerperants data, and the Beetles data. It should be noted that the covariance matrix of and is based on large sample theory. Hence, is replaced by when calculating the CIs for all the methods.

4.2 Real-World Datasets

The Hewlett data are a classical data example with nine dose levels, a relatively moderate sample size, and a steep response curve. This dataset was first used in a paper by Abdelbasit and Plackett (1983) where the authors stated that it was obtained from a personal communication with Hewlett without any more details. Subsequently, it was further studied in several related papers (Sitter and Wu, 1993; Faraggi et al., 2003; Paige et al., 2011). The Puerperants data were collected to estimate the of levobupivacaine for labor analgesia (Sia et al., 2005). This dataset has five dose levels, and a steep response curve, but a relatively small sample size. In these data, 50 parturients in early labor were randomly assigned to receive one of the five doses studied. Effective analgesia is defined as a pain score (0-100 visual analog scale) of less than 10 within 15 min of injection, lasting for 45 min or more after the induction of analgesia. The Beetles data studied the toxicity of beetles to insecticide. This dataset has six dose levels, and a moderate sample size, but a relatively flat response curve. It is taken from Hewlett and Plackett (1950), and has been further discussed by Zelterman (1999) and Faraggi et al. (2003). In this study, groups of beetles were exposed to various concentrations of the insecticide. After a short period, the number of deaths within each group was reported. Web Tables 18-20, respectively, provide the details of all three datasets. Further, we fit the standard logistic regression to each of the three datasets, and the summary statistics are displayed in Table 3.

4.3 Further Simulation Study Based on the Hewlett Data

For the real data described in Section 4.2, both and are estimated from binary data, but such case has not been studied in the simulation study in Section 3. It is known that the mean and the variance of the logistic regression coefficient are correlated especially in small samples. Therefore, based on the Hewlett data, we further conduct an extensive simulation study, with various slopes and sample sizes, to investigate the performance of the PF method in this scenario, and the details are provided in Web Appendix  C. Based on the new simulation results in Web Tables 21 and 22, we conclude that the PF method has the overall best performance in terms of coverage, width, and location, and thus is still recommended in the binary data case (see description from Web Appendix  C).

4.4 Data Analysis

Table 4 displays the CIs estimated from all four methods based on the Hewlett data, the Puerperants data, and the Beetles data. In addition, we also evaluate the performance of all the methods in the scenario where the parameters ⁠, ⁠, γ, and β are taken or estimated from each of the datasets. From Table 4, we observe that both the PF and Fieller methods generally have a valid CP, but the former can provide a narrower CW than the latter, when the sample size is small or the response curve is flat. Meanwhile, it is not surprising to see that the delta method has the narrowest CW among all the methods. However, this method should be used with caution, because it may underestimate the CP even in moderate samples. For the Beetles data, which has a relatively flat response curve, although the delta method seems to have the best performance, this method can seriously overestimate the CP and may not be reliable as discussed in Web Appendix  C. Further, it can be seen that the PB method has the worst performance and thus is also not recommended for practical use.

In conclusion, the PF method performs well as long as the sample size at each dose is not very small and the response curve is neither extremely flat nor extremely steep. However, if such conditions are violated, all four methods may not work well. This is probably because the covariance matrix for the numerator and the denominator cannot be estimated accurately. Under this circumstance, the likelihood-based CI may be more reliable (Williams, 1986). It is therefore of interest to generalize the PF method to the penalized likelihood method. We leave this for future research.

5 Discussion

In this article, we have developed the penalized Fieller method to construct the CI for the ratio estimate. Like the Fieller method, the proposed approach is based on the inverting of the t statistic, and thus it is an analytic approach and naturally allows an asymmetric CI. Further, to overcome the unbounded issue of the Fieller CI, we adopt a penalized likelihood approach to estimate the denominator while adjusting the numerator accordingly to reduce the bias of the t statistic. By using a suitable penalty parameter, the PF method always produces a bounded CI. Moreover, we show theoretically that the PF method is second-order correct under the bivariate normality assumption, better than other bounded methods, which are all first-order correct. Simulation results demonstrate that the PF method has good performance in terms of controlling the CP and the CW. Hence, we recommend using the PF method for its robustness and analytic advantage.

Note that we obtain the penalized estimator by directly plugging into it. By using this strategy, we can get a much simpler expression of ⁠, and therefore simplify the subsequent statistical inferences as well. We can also obtain the joint estimates and by maximizing the full penalized likelihood, and the details are provided in Web Appendix  D. It is seen that is much more complicated but will converge to as ⁠. Note that the forms of and only differ at the variance estimate. For the variance estimate ⁠, it is related to the penalty parameter λ only at the third order. This fact indicates that it makes little difference to consider a penalty when estimating v₂. Therefore, we directly plug into for mathematical convenience.

It should be emphasized that the PF method is not an exact method. That is, if the sample size is extremely small, and the denominator has a very large variation (the case in which the Fieller method has a very large UP), then the PF method still has the possibility of underestimating the CP as shown in the simulation results. Although the Fieller method is an exact method in the bivariate normality case, and is generally guaranteed to achieve the nominal level, this is at the potential cost of having intervals of infinite length. Further, the Fieller method ignores the information that the denominator cannot be 0. By excluding the case of ⁠, it is not surprising that the PF method has a narrower CW than the Fieller method. Notably, the denominator being not statistically different from 0 also means that the CI of the denominator covers 0. Because 0 is a singular point for the denominator of a ratio, it is expected that the Fieller CI will be unbounded in this case. However, the true value of the ratio must be a real number. Therefore, the unbounded interval is not as informative as the bounded interval because the latter always provides finite length and is much more common in practice. On the other hand, the PF method obtains the bounded CI by shrinking the CI of the denominator away from 0 through penalty, while adjusting the CI of the numerator simultaneously.

Note that the PF method is second-order correct only in the bivariate normality case. Otherwise, all the existing methods, including the PF and Fieller methods, are first-order correct. Like the Fieller method, despite the fact that the PF method is proposed based on the bivariate normality assumption, we have demonstrated its usefulness in a wide range of scenarios through simulation study and real data application. These scenarios include the ratio of means from both the bivariate normal and nonnormal distributions as well as the ratio of the nonlinear regression coefficients. However, in certain cases, when the underlying distribution is very skewed, the PF method may not work well, just like almost all the existing methods. This is probably due to two reasons. First, the mean estimator may need a much larger sample size in order to converge to a normal distribution. Second, using the mean to estimate the location parameter is no longer appropriate. In the case of very skewed distributions, the CI based on inverting of the nonparametric statistic may be more useful, but this needs further research. On the other hand, it is known that the MOVER-R (method of variance estimate recovery for a ratio) is a generalization of the Fieller method to the nonnormal case (Donner and Zou, 2012; Newcombe, 2016). So, developing the penalized MOVER-R is another possible direction for future research.

Recall the fact that we only use the penalized estimate for the denominator by shrinking it away from 0. Because the numerator generally has the possibility to be 0, we do not penalize both the denominator and the numerator simultaneously. For example, considering the median response dose ⁠, it is possible that the numerator β₀ can equal 0. This corresponds to the scenario of ⁠. Because γ is generally on a logarithmic scale, the original median response dose will be 1. In some cases, if we know the information that the numerator cannot be 0, we can use the penalized estimate for both the numerator and denominator. By using more information, it is possible to get a narrower CW. However, how to obtain a second-order correct CI in this situation is not at all clear.

Finally, we use the penalized ratio estimate to estimate r in ⁠. Note that there exists other point estimates for r that we can select to plug into (Tin, 1965), which might be better than ⁠. However, further studies are needed to investigate their performance, especially when the denominator does not have adequate power to reject being 0.

Data Availability Statement

The data that support the findings in this paper are available in the Supporting Information of this article.

Acknowledgments

The authors would like to thank the associate editor, the two anonymous reviewers, and the co-editor for their constructive comments which greatly improved the presentation of this article. The authors also appreciate the help of Raymond J. Carroll, PhD for his critical reading of the original version of the manuscript. This work is partially supported by National Institutes of Health grants R03-DE024198 and R03-DE025646.

References

Appendix A: Maximum Point of the Penalized Log-Likelihood

Treating v₂ as known and taking the first derivative with respect to μ₂, we obtain

where and ⁠. It is easy to see that regardless of the sign of ⁠. Therefore, the penalized log-likelihood will increase when ⁠. That is, the penalized log-likelihood can only attain its maximum value at the point x₁ or x₂. Further,

If ⁠, we have ⁠, which implies ⁠. However, if ⁠, we obtain ⁠. Hence, the penalized log-likelihood attains its maximum value at the point ⁠.

APPENDIX B: Confidence Width of the penalized Fieller Method When λ → + ∞

Let denote the CW of the PF CI. According to Section 2.6, it is seen that when λ is greater than or equal to ⁠, where ⁠, ⁠, and ⁠. In order to avoid confusion, we use the notation to represent the order with respect to λ. It is easy to see that ⁠, ⁠, and ⁠. Hence, ⁠. That is, when ⁠, we have ⁠, and the rate of convergence is of order ⁠.

Appendix C: Proof of Theorem 1

Let and denote the UP of the PF and Fieller methods, respectively. According to Fieller's theorem, the PF CI is unbounded if and only if

Substituting into the above inequality, we have

Therefore,

Specially, can always be guaranteed if and only if the penalty parameter ⁠. This completes the proof.

Appendix D: Proof of Theorem 2

Under the bivariate normality assumption, the Fieller method is known as the only exact method (Koschat, 1981). Without loss of generality, let and ⁠, respectively, denote the confidence limit for the PF and Fieller methods that satisfy

When ⁠, both and are lower limits; otherwise, they are upper limits. To prove the theorem, we need to show that ⁠.

Because μ₂ is a nonzero constant, we have and ⁠. From Equation (A1), we know that ⁠. Therefore,

Using the above equation, we rewrite the Equation (A1) as

Note that

The first part of the above equation is and the second part can be expressed as ⁠, where

Equation (A2) minus Equation (A3) becomes

By the fact ⁠, we can see ⁠, which completes the proof.