The Owen T-function
The OwenQ package also provides an implementation of the Owen T-function, under the name OwenT. This is a port of the function owens_t of boost, the peer-reviewed collection of C++ libraries.
Let Z∼N(0,1) and X∼χ2ν be two independent random variables. For real numbers δ1 and δ2, define the two random variables T1=Z+δ1√X/νandT2=Z+δ2√X/ν.
Both T1 and T2 follow a non-central Student distribution. The number of degrees of freedom is ν for each of them, and their respective non-centrality parameters are δ1 and δ2 respectively.
Owen (1965) studied the distribution of the pair (T1,T2).
The four Owen cumulative functions are O1(ν,t1,t2,δ1,δ2)=Pr
Owen provided an efficient way to evaluate these functions when \nu is an integer number. Owen’s algorithms are implemented in the OwenQ package.
For \delta_1 > \delta_2, these four functions are implemented in the OwenQ package under the respective names powen1, powen2, powen3 and powen4. For general values of \delta_1 and \delta_2, they are implemented under the respective names psbt1, psbt2, psbt3 and psbt4.
Owen (1965) also provided an algorithm to evaluate the cumulative distribution function of a univariate non-central Student distribution with an integer number of degrees of freedom. This evaluation is performed by the function ptOwen of the OwenQ package.
It is known that the pt function is not reliable when the non-centrality parameter ncp is large. Below we compare the values given by ptOwen and pt to the value given by Wolfram|Alpha (returned by the command N[CDF[NoncentralStudentTDistribution[4,70],80]]):
When q=delta, the value of ptOwen(q, nu, delta) should go to 0.5 as nu increases to infinity. The examples below show the failure of this expectation when nu is too large.
ptOwen(q=50, nu=3500, delta=50)
## [1] 0.4986697
ptOwen(q=50, nu=3600, delta=50)
## [1] 0.4989289
ptOwen(q=50, nu=3650, delta=50)
## [1] 0.4522949
ptOwen(q=50, nu=3660, delta=50)
## [1] 0.3349762
ptOwen(q=50, nu=3670, delta=50)
## [1] 0.7607795
ptOwen(q=50, nu=3680, delta=50)
## [1] 0Since all Owen’s algorithms are somehow similar to the algorithm evaluating ptOwen, we can suspect that the other ones also suffer from certain limitations. In the other vignette, we investigate the reliability and the limitations of OwenQ.
In order to do some comparisons, and thanks to the BH package, we have ported the boost implementation of the cumulative Student distribution to OwenQ. It is evaluated by the internal function pt_boost. We concluded that this function is highly reliable. In particular it does not suffer from the failures of ptOwen we have just shown:
OwenQ:::pt_boost(q=50, nu=3500, delta=50)
## [1] 0.4986697
OwenQ:::pt_boost(q=50, nu=3600, delta=50)
## [1] 0.4987041
OwenQ:::pt_boost(q=50, nu=3650, delta=50)
## [1] 0.4987206
OwenQ:::pt_boost(q=50, nu=3660, delta=50)
## [1] 0.4987238
OwenQ:::pt_boost(q=50, nu=3670, delta=50)
## [1] 0.4987271
OwenQ:::pt_boost(q=50, nu=3680, delta=50)
## [1] 0.4987303The Owen distribution intervenes in the calculation of the power of equivalence tests.
Assume a statistical model given by a sample y_i \sim_{\text{iid}} {\cal N}(\mu, \sigma^2) for i=1, \ldots, n. We want to demonstrate that, up to a given confidence level, the mean \mu belongs to a certain interval [\Delta_1, \Delta_2]. In other words, we are interested in the alternative hypothesis H_1\colon\{\Delta_1 \leq \mu \leq \Delta_2\}.
Consider the usual 100(1-2\alpha)\%-confidence interval about \mu: \left[\bar y - t^\ast_{n-1}(\alpha)\frac{\hat\sigma}{\sqrt{n}}, \, \bar y + t^\ast_{n-1}(\alpha)\frac{\hat\sigma}{\sqrt{n}} \right].
The H_1 hypothesis is accepted at level \alpha when this interval falls into the interval [\Delta_1, \Delta_2].
This can be written as follows: T_1 := \frac{\bar y - \Delta_1}{\hat\sigma/\sqrt{n}} \geq t^\ast_{n-1}(\alpha) \quad \text{and} \quad T_2 := \frac{\bar y - \Delta_2}{\hat\sigma/\sqrt{n}} \leq - t^\ast_{n-1}(\alpha).
Observe that T_1 = \frac{z - \delta_1}{\dfrac{\sqrt{n-1}\hat\sigma/\sigma}{\sqrt{n-1}}} where z = \frac{\sqrt{n}}{\sigma}(\bar y - \mu) \sim {\cal N}(0,1) \quad \text{and} \quad \delta_1 = \frac{\mu - \Delta_1}{\sigma/\sqrt{n}}.
By reasoning in the same way for T_2, we find that the pair (T_1, T_2) follows the Owen distribution with degrees of freedom \nu = n-1, and non-centrality parameters \delta_1 given above and \delta_2 = \tfrac{\mu - \Delta_2}{\sigma/\sqrt{n}}.
Therefore the power of the test - i.e. the probability to accept H_1 - is given by
O_4\bigl(n-1, t^\ast_{n-1}(\alpha), -t^\ast_{n-1}(\alpha), \delta_1, \delta_2\bigr),
and then can be evaluated thanks to powen4.
The result of the equivalence test is said to be inconclusive when only one of the bounds of the confidence interval falls into the interval [\Delta_1, \Delta_2].
The probability to get an inconclusive result can be obtained with the OwenQ package. We show and check that below with the help of simulations.
Delta1 <- -2; Delta2 <- 2
mu <- 1; sigma <- 6; n <- 30L
alpha <- 0.05
nsims <- 1e6L
equivalence <- inconclusive <- numeric(nsims)
for (i in 1L:nsims) {
y <- rnorm(n, mu, sigma)
CI <- t.test(x = y, conf.level = 1-2*alpha)$conf.int
equivalence[i] <- (CI[1] > Delta1) && (CI[2] < Delta2)
inconclusive[i] <- ((CI[1] < Delta1) && (CI[2] > Delta1)) ||
((CI[1] < Delta2) && (CI[2] > Delta2))
}dof <- n-1
q <- qt(1-alpha, dof)
se <- sqrt(1/n)*sigma
delta1 <- (mu-Delta1)/se; delta2 <- (mu-Delta2)/se
# probability to get equivalence
mean(equivalence)
## [1] 0.092927
powen4(dof, q, -q, delta1, delta2)
## [1] 0.09300963
# probability to get inconclusive
mean(inconclusive)
## [1] 0.901587
ptOwen(q, dof, delta2) - ptOwen(-q, dof, delta1) - powen4(dof, q, -q, delta1, delta2)
## [1] 0.90139Now consider two independent samples x_i \sim_{\text{iid}} {\cal N}(\mu, \sigma^2) for i=1, \ldots, n_1. and y_i \sim_{\text{iid}} {\cal N}(\nu, \sigma^2) for i=1, \ldots, n_2 and the alternative hypothesis H_1\colon\bigl\{|\mu-\nu| \leq \Delta\bigr\}. The power of this test is evaluated by the function powerTOST below. The parameter delta0 corresponds to the difference \mu-\nu.
The OwenQ package also provides an implementation of the Owen T-function, under the name OwenT. This is a port of the function owens_t of boost, the peer-reviewed collection of C++ libraries.
The Owen cumulative functions are based on the two Owen Q-functions Q_1(\nu, t, \delta, R) = \frac{1}{\Gamma\left(\frac{\nu}{2}\right)2^{\frac12(\nu-2)}} \int_0^R \Phi\left(\frac{tx}{\sqrt{\nu}}-\delta\right) x^{\nu-1} e^{-\frac{x^2}{2}} \mathrm{d}x and Q_2(\nu, t, \delta, R) = \frac{1}{\Gamma\left(\frac{\nu}{2}\right)2^{\frac12(\nu-2)}} \int_R^\infty \Phi\left(\frac{tx}{\sqrt{\nu}}-\delta\right) x^{\nu-1} e^{-\frac{x^2}{2}} \mathrm{d}x.
They are implemented in the OwenQ package under the respective names OwenQ1 and OwenQ2 (for integer values of \nu).
Consider the statistical model given by a sample y_i \sim_{\text{iid}} {\cal N}(\mu, \sigma^2) for i=1, \ldots, n.
An equal-tailed (p,\alpha)-tolerance interval is an interval containing at least 100p\% of the “center data” with 100(1-\alpha)\% confidence.
The natural choice for such an interval is, as shown by Owen (1965), \bigl[\bar y - k_e \hat\sigma, \bar y + k_e \hat\sigma] where k_e satisfies O_2(n-1, k_e\sqrt{n}, -k_e\sqrt{n}, \delta, -\delta) = 1-\alpha where \delta = \sqrt{n}z_{(1+p)/2}.
Therefore, the tolerance factor k_e can be determined with the help of the powen2 function of the OwenQ package. But it is more efficient to use the function spowen2 for this purpose; spowen2(nu, t, delta) has the same value as powen2(nu, t, -t, delta, -delta) but it is evaluated more efficiently.
p <- 0.9; alpha <- 0.05
n <- 100
delta <- sqrt(n)*qnorm((1+p)/2)
uniroot(function(ke) spowen2(n-1, ke*sqrt(n), delta) - (1-alpha),
lower=qnorm(1-alpha), upper=5, extendInt = "upX", tol=1e-9)$root
## [1] 1.981513The k_e factors are tabulated in Krishnamoorthy & Mathew’s book.
The OwenQ package provides three internal functions, iOwenQ1, iOwenQ2 and ipowen4, which respectively perform the evaluation of Q_1, Q_2 and O_4 by numerical integration using the RcppNumerical package. They can be called with a positive non-integer value of \nu. The evaluation of O_4 with ipowen4 is correct only for t_1 > t_2 and \delta_1 > \delta_2.
Owen, D. B. (1965). A special case of a bivariate noncentral t-distribution. Biometrika 52, 437-446.
Krishnamoorthy & Mathew (2009). Statistical Tolerance Regions. Wiley.