Purpose

The purpose of this vignette is to assess the correctness of the functions of the OwenQ package.

As said in the main vignette, the fourth Owen cumulative function O₄, implemented under the name powen4 allows to evaluate the power of equivalence tests. This is perhaps the main application of the OwenQ package and we will particularly focus on this situation.

Our strategy runs as follows:

• To test the OwenT function, we will compare some values it returns to those returned by other implementations. Note that the other functions of the OwenQ package rely on OwenT when the number of degrees of freedom is odd. Therefore, by testing these other functions for an odd number of degrees of freedom, we implicitely test the OwenT function.

• We will compare a couple of values returned by OwenQ1 and OwenQ2 to the values obtained by numerical integration in Wolfram|Alpha.

• We will test powen4 by performing some power calculations relying on this function and comparing the results with the ones provided by SAS. We will also perform these power calculations with the help of the ipowen4 function, which evaluates the O₄ function by a powerful numerical integration, and we will compare the results.

• These power calculations will be performed for 100 different scenarios, and we will conclude that powen4 is reliable for these scenarios. Then, to test the other functions of the OwenQ package, we will check that some mathematical relations hold for these 100 scenarios.

The folder tests/testthat of the OwenQ package contains many comparisons with Wolfram|Alpha. The results we obtain with OwenQ are always the same as the ones obtained with Wolfram|Alpha up to a tolerance of at least 10 decimal digits.

Owen T-function (OwenT)

The Owen T-function is implemented under the name OwenT. This is a port of the function owens_t of the C++ special_functions library included in the boost set of libraries. Some details about the C++ implementation are available here.

The libraries provided by boost are peer-reviewed, and this is enough to trust the reliability of OwenT. Nevertheless, below we compare the results of OwenT to the six values given in Patefield’s article, up to 14 decimal digits. We have also checked that Wolfram|Alpha gives these values.

testData <- data.frame(
  h = c(0.0625, 6.5, 7, 4.78125, 2, 1), 
  a = c(0.25, 0.4375, 0.96875, 0.0625, 0.5, 0.9999975),
  Patefield = c(3.8911930234701e-02, 2.0005773048508e-11, 6.3990627193899e-13, 1.0632974804687e-07, 8.6250779855215e-03, 6.6741808978229e-02),
  OwenT = numeric(6)
)
for(i in 1:nrow(testData)){
  testData$OwenT[i] <- OwenT(testData$h[i], testData$a[i])
}
print(testData, digits=14)
##         h         a           Patefield               OwenT
## 1 0.06250 0.2500000 3.8911930234701e-02 3.8911930234701e-02
## 2 6.50000 0.4375000 2.0005773048508e-11 2.0005773048508e-11
## 3 7.00000 0.9687500 6.3990627193899e-13 6.3990627193899e-13
## 4 4.78125 0.0625000 1.0632974804687e-07 1.0632974804687e-07
## 5 2.00000 0.5000000 8.6250779855215e-03 8.6250779855215e-03
## 6 1.00000 0.9999975 6.6741808978229e-02 6.6741808978229e-02

We get the same results with OwenT.

Owen Q-functions: comparing a couple of values

The two Owen Q-functions Q₁ and Q₂ are defined by: $$ 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 ν, following Owen’s algorithm (1965).

In Wolfram|Alpha, these functions are not available, but we can evaluate them by numerical integration.

Below we compare two values of OwenQ1 to the ones returned by Wolfram|Alpha.

• ν = 3, t = 3, δ = 2, R = 5 (link to Wolfram)

# wolfram: Integrate[(1+Erf[(3*x/Sqrt[3]-2)/Sqrt[2]])*x^(3-1)*Exp[-x^2/2],{x,0,5}]/2/Gamma[3/2]/2^((3-2)/2)
OwenQ1(3, 3, 2, 5)
## [1] 0.6800117

Our value rounded to 6 digits is the same as the one given by Wolfram.

• ν = 1000, t = 3, δ = 2, R = 30 (link to Wolfram)

# wolfram: Integrate[(1+Erf[(3*x/Sqrt[1000]-2)/Sqrt[2]])*x^(1000-1)*Exp[-x^2/2],{x,0,30}]/2/Gamma[1000/2]/2^((1000-2)/2)
print(OwenQ1(1000, 3, 2, 30), digits=16)
## [1] 0.008518809463589428

The two values are the same up to 13 digits.

Now we compare two values of OwenQ2 to the ones returned by Wolfram|Alpha.

• ν = 3, t = 3, δ = 2, R = 5 (link to Wolfram)

# wolfram: Integrate[(1+Erf[(3*x/Sqrt[3]-2)/Sqrt[2]])*x^(3-1)*Exp[-x^2/2],{x,5,Infinity}]/2/Gamma[3/2]/2^((3-2)/2)
OwenQ2(3, 3, 2, 5)
## [1] 1.54405e-05

The two values are identical.

• ν = 1000, t = 3, δ = 2, R = 5 (link To Wolfram)

# wolfram: Integrate[(1+Erf[(3*x/Sqrt[1000]-2)/Sqrt[2]])*x^(1000-1)*Exp[-x^2/2],{x,5,Infinity}]/2/Gamma[1000/2]/2^((1000-2)/2)
print(OwenQ2(1000, 3, 2, 5), digits=16)
## [1] 0.8406201459600969

The two values are identical up to 13 digits.

Fourth Owen cumulative function O₄ (powen4)

As seen in the other vignette, the fourth Owen cumulative function O₄, implemented under the name powen4, can be used to evaluate the power of equivalence tests.

The powerTOST function below returns the power of the equivalence test for a so-called parallel design with equal variances, when considering the alternative hypothesis H₁: {−Δ < δ₀ < Δ}, where δ₀ denotes the difference between the two means. This function takes as arguments the significance level α, the difference δ₀ between the two means, the threshold Δ, the common standard deviation σ of the two samples, and the two sample sizes n₁ and n₂.

powerTOST <- function(alpha, delta0, Delta, sigma, n1, n2, algo=2) {
  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, algo=algo)
}

As you can see, the powen4 function has an argument algo. This is also the case for the other Owen functions, except for ptOwen. The algo argument can take two values, 1 or 2. Th default value is algo=2, and as we will see later, the evaluation is more reliable with this option. The value of algo corresponds to a small difference in the algorithm. With algo=1, the evaluation is a bit faster. And we will see that powerTOST
is reliable for algo=1 when the value of n1+n2 is not too large.

power <- numeric(nrow(SAS))
for (i in 1:nrow(SAS)) {
  power[i] <-
    powerTOST(
      alpha = SAS$alpha[i],
      delta0 = SAS$delta0[i],
      Delta = SAS$Delta[i],
      sigma = SAS$sigma[i],
      n1 = SAS$n1[i],
      n2 = SAS$n2[i]
    )
}

identical(round(power,5), SAS$powerSAS)
## [1] TRUE

power <- numeric(nrow(SAS))
for (i in 1:nrow(SAS)) {
  power[i] <-
    powerTOST(
      alpha = SAS$alpha[i],
      delta0 = SAS$delta0[i],
      Delta = SAS$Delta[i],
      sigma = SAS$sigma[i],
      n1 = SAS$n1[i],
      n2 = SAS$n2[i],
      algo = 1
    )
}
identical(round(power,5), SAS$powerSAS)
## [1] TRUE

ipowerTOST <- 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)
  OwenQ:::ipowen4(dof, q, -q, delta1, delta2)
}

alpha <- 0.05; delta0 <- 0; Delta <- 5
sigma <- 110
n1 <- n2 <- 2500
powerTOST(alpha, delta0, Delta, sigma, n1, n2)
## [1] 4.523596e-05
ipowerTOST(alpha, delta0, Delta, sigma, n1, n2)
## [1] 4.523596e-05
## attr(,"err_est")
## [1] 5.022201e-19
## attr(,"err_code")
## [1] 0

sigma <- 152
n1 <- n2 <- 5000
powerTOST(alpha, delta0, Delta, sigma, n1, n2)
## [1] 0.003612374
ipowerTOST(alpha, delta0, Delta, sigma, n1, n2)
## [1] 0.003612374
## attr(,"err_est")
## [1] 4.010541e-17
## attr(,"err_code")
## [1] 0

power <- ipower <- numeric(nrow(SAS))
for (i in 1:nrow(SAS)) {
  power[i] <-
    powerTOST(
      alpha = SAS$alpha[i],
      delta0 = SAS$delta0[i],
      Delta = SAS$Delta[i],
      sigma = SAS$sigma[i],
      n1 = SAS$n1[i],
      n2 = SAS$n2[i]
    )
  ipower[i] <-
    ipowerTOST(
      alpha = SAS$alpha[i],
      delta0 = SAS$delta0[i],
      Delta = SAS$Delta[i],
      sigma = SAS$sigma[i],
      n1 = SAS$n1[i],
      n2 = SAS$n2[i]
    )
}
identical(round(power, 10), round(ipower, 10))
## [1] TRUE

power <- numeric(nrow(SAS))
for (i in 1:nrow(SAS)) {
  power[i] <-
    powerTOST(
      alpha = SAS$alpha[i],
      delta0 = SAS$delta0[i],
      Delta = SAS$Delta[i],
      sigma = SAS$sigma[i],
      n1 = SAS$n1[i],
      n2 = SAS$n2[i], 
      algo = 1
    )
}
identical(round(power, 10), round(ipower, 10))
## [1] TRUE

Other Owen cumulative functions

We previously validated the powen4 function for the 100 scenarios of the dataset SAS.

Now we test the functions powen1, powen2 and powen3 by checking the equality O₁(ν, t₁, t₂, δ₁, δ₂) + O₂(ν, t₁, t₂, δ₁, δ₂) + O₃(ν, t₁, t₂, δ₁, δ₂) + O₄(ν, t₁, t₂, δ₁, δ₂) = 1.

We restrict our attention to default setting algo=2. We check the above equality for the 100 scenarios of the dataset SAS.

f <- 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)
  powen1(dof, q,-q, delta1, delta2) + powen2(dof, q,-q, delta1, delta2) + 
    powen3(dof, q,-q, delta1, delta2) + powen4(dof, q,-q, delta1, delta2)
}
test <- numeric(nrow(SAS))
for (i in 1:nrow(SAS)) {
  test[i] <-
    f(
      alpha = SAS$alpha[i],
      delta0 = SAS$delta0[i],
      Delta = SAS$Delta[i],
      sigma = SAS$sigma[i],
      n1 = SAS$n1[i],
      n2 = SAS$n2[i]
    )
}
all(abs(test-1) < 1e-14)
## [1] TRUE

We find that each of the 100 sums is equal to 1 up to a tolerance of 14 digits.

First Owen Q-function Q₁ (OwenQ1)

We previously validated the powen4 function (implementation of the fourth Owen cumulative function O₄) for the 100 scenarios of the dataset SAS.
Now, to test OwenQ1, we will use the following equality:

$$ O_4(\nu, t_1, t_2, \delta_1, \delta_2) = Q_1\left(\nu, t_2, \delta_2, \frac{\sqrt{\nu}(\delta_1-\delta_2)}{t_1-t_2}\right) - Q_1\left(\nu, t_1, \delta_1, \frac{\sqrt{\nu}(\delta_1-\delta_2)}{t_1-t_2}\right). $$

We use this formula to perform the power calculations with OwenQ1 instead of powen4.

powerTOST2 <- 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)
  R <- sqrt(dof)*(delta1 - delta2)/q/2
  OwenQ1(dof, -q, delta2, R) - OwenQ1(dof, q, delta1, R)
}
power2 <- numeric(nrow(SAS))
for (i in 1:nrow(SAS)) {
  power2[i] <-
    powerTOST2(
      alpha = SAS$alpha[i],
      delta0 = SAS$delta0[i],
      Delta = SAS$Delta[i],
      sigma = SAS$sigma[i],
      n1 = SAS$n1[i],
      n2 = SAS$n2[i]
    )
}
all(abs(power - power2) < 1e-9)
## [1] TRUE

The values rounded to 9 digits are the same as the ones we previously obtained with powen4.

Second Owen Q-function Q₂ (OwenQ2)

We previously validated powen2 for the 100 scenarios of the dataset SAS. Now we test the OwenQ2 function by checking the equality $$ O_2(\nu, t_1, t_2, \delta_1, \delta_2) = Q_2\left(\nu, t_1, \delta_1, \frac{\sqrt{\nu}(\delta_1-\delta_2)}{t_1-t_2}\right) - Q_2\left(\nu, t_2, \delta_2, \frac{\sqrt{\nu}(\delta_1-\delta_2)}{t_1-t_2}\right). $$

We check this equality for the 100 scenarios of the dataset SAS.

g <- 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)
  R <- sqrt(dof)*(delta1 - delta2)/q/2
  x <- OwenQ2(dof, q, delta1, R) - OwenQ2(dof, -q, delta2, R)
  y <- powen2(dof, q, -q, delta1, delta2)
  x - y
}
test <- numeric(nrow(SAS))
for (i in 1:nrow(SAS)) {
  test[i] <-
    g(
      alpha = SAS$alpha[i],
      delta0 = SAS$delta0[i],
      Delta = SAS$Delta[i],
      sigma = SAS$sigma[i],
      n1 = SAS$n1[i],
      n2 = SAS$n2[i]
    )
}
all(abs(test) < 1e-15)
## [1] TRUE

We find that each of the 100 equalities hold true up to a tolerance of 15 digits.

Student cumulative distribution function (ptOwen)

We finally test ptOwen, the implementation of the cumulative distribution function of the noncentral Student distribution.
If T₁ follows the noncentral Student distribution with number of degrees of freedom ν and noncentrality parameter δ₁, then for any t₁ the equality Pr (T₁ ≤ t₁) = O₁(ν, t₁, t₂, δ₁, δ₂) + O₂(ν, t₁, t₂, δ₁, δ₂) holds for any t₂ and δ₂. We check this equality for the 100 scenarios of the dataset SAS.

h <- 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)
  x <- ptOwen(q, dof, delta1)
  y <- powen1(dof, q, -q, delta1, delta2) + powen2(dof, q, -q, delta1, delta2)
  x - y
}
test <- numeric(nrow(SAS))
for (i in 1:nrow(SAS)) {
  test[i] <-
    h(
      alpha = SAS$alpha[i],
      delta0 = SAS$delta0[i],
      Delta = SAS$Delta[i],
      sigma = SAS$sigma[i],
      n1 = SAS$n1[i],
      n2 = SAS$n2[i]
    )
}
all(abs(test) < 1e-15)
## [1] TRUE

For each scenario, the equality holds up to a tolerance of 15 digits.

References

• Patefield, M. (2000). Fast and Accurate Calculation of Owen’s T Function. Journal of Statistical Software 5 (5), 1-25.

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

• Wolfram Alpha LLC. 2017. Wolfram|Alpha. (access July 17, 2017).