In stla/OwenQ: Owen Q-Function

library(OwenQ)
knitr::opts_chunk$set(collapse=TRUE)

The Owen distribution

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.

Non-central Student distribution

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)
pt(q=1, df=3, ncp=2)

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]]):

p1 <- pt(q=80, df=4, ncp=70)
p2 <- ptOwen(q=80, nu=4, delta=70)
wolfram <- 0.54742763380700947685
p1 - wolfram
p2 - wolfram

Limitations

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)
ptOwen(q=50, nu=3600, delta=50)
ptOwen(q=50, nu=3650, delta=50)
ptOwen(q=50, nu=3660, delta=50)
ptOwen(q=50, nu=3670, delta=50)
ptOwen(q=50, nu=3680, delta=50)

Since 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)
OwenQ:::pt_boost(q=50, nu=3600, delta=50)
OwenQ:::pt_boost(q=50, nu=3650, delta=50)
OwenQ:::pt_boost(q=50, nu=3660, delta=50)
OwenQ:::pt_boost(q=50, nu=3670, delta=50)
OwenQ:::pt_boost(q=50, nu=3680, delta=50)

Application to equivalence testing

The Owen distribution intervenes in the calculation of the power of equivalence tests.

One sample

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.

Inconclusive equivalence test

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)
powen4(dof, q, -q, delta1, delta2)
# probability to get inconclusive
mean(inconclusive)
ptOwen(q, dof, delta2) - ptOwen(-q, dof, delta1) - powen4(dof, q, -q, delta1, delta2)

Parallel design

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

h <- seq(-3, 3, length.out=30)
a <- seq(-5, 5, length.out=30)
z <- outer(h, a, Vectorize(OwenT))
oldpar <- par(mar=c(0,2,0,0))
persp(h, a, z, theta=30, phi=30, expand=0.5, zlab="Owen T", ticktype = "detailed", col = "lightblue")
par(oldpar)

The Owen $Q$-functions

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

Application: Equal-tailed tolerance interval

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

The $k_e$ factors are tabulated in Krishnamoorthy & Mathew's book.

Internal functions: numerical integration

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

References

• Owen, D. B. (1965). A special case of a bivariate noncentral t-distribution. Biometrika 52, 437-446.

• Krishnamoorthy & Mathew (2009). Statistical Tolerance Regions. Wiley.

stla/OwenQ documentation built on April 24, 2023, 9:49 p.m.