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 \sim {\cal N}(0,1)\) and \(X \sim \chi^2_\nu\) be two independent random variables. For real numbers \(\delta_1\) and \(\delta_2\), define the two random variables \[ T_1 = \frac{Z+\delta_1}{\sqrt{X/\nu}} \quad \text{and} \quad\; T_2 = \frac{Z+\delta_2}{\sqrt{X/\nu}}. \]
Both \(T_1\) and \(T_2\) follow a non-central Student distribution. The number of degrees of freedom is \(\nu\) for each of them, and their respective non-centrality parameters are \(\delta_1\) and \(\delta_2\) respectively.
Owen (1965) studied the distribution of the pair \((T_1, T_2)\).
The four Owen cumulative functions are \[ \begin{align} O_1(\nu, t_1, t_2, \delta_1, \delta_2) & = \Pr(T_1 \leq t_1, T_2 \leq t_2), \\ O_2(\nu, t_1, t_2, \delta_1, \delta_2) & = \Pr(T_1 \leq t_1, T_2 \geq t_2), \\ O_3(\nu, t_1, t_2, \delta_1, \delta_2) & = \Pr(T_1 \geq t_1, T_2 \geq t_2), \\ O_4(\nu, t_1, t_2, \delta_1, \delta_2) & = \Pr(T_1 \geq t_1, T_2 \leq t_2). \end{align} \]
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.
ptOwen(q=1, nu=3, delta=2)
## [1] 0.1573494
pt(q=1, df=3, ncp=2)
## [1] 0.1573494It 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]]):
p1 <- pt(q=80, df=4, ncp=70)
p2 <- ptOwen(q=80, nu=4, delta=70)
wolfram <- 0.54742763380700947685
p1 - wolfram
## [1] 0.02268711
p2 - wolfram
## [1] 3.330669e-16When 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.092532
powen4(dof, q, -q, delta1, delta2)
## [1] 0.09300963
# probability to get inconclusive
mean(inconclusive)
## [1] 0.901865
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\).
powerTOST <- function(alpha, delta0, Delta, sigma, n1, n2) {
se <- sqrt(1/n1 + 1/n2) * sigma
delta1 <- (delta0 + Delta) / se
delta2 <- (delta0 - Delta) / se
dof <- n1 + n2 - 2
q <- qt(1 - alpha, dof)
powen4(dof, q, -q, delta1, delta2)
}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.