## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set(collapse=TRUE) library(OwenQ) ## ----------------------------------------------------------------------------- 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) ## ----------------------------------------------------------------------------- # 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) ## ----------------------------------------------------------------------------- # 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) ## ----------------------------------------------------------------------------- # 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) ## ----------------------------------------------------------------------------- # 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) ## ----------------------------------------------------------------------------- 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) } ## ----echo=FALSE--------------------------------------------------------------- SAS <- structure(list(alpha = c(0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.01, 0.01, 0.01, 0.01, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.05, 0.05, 0.01, 0.01, 0.05, 0.05), delta0 = c(0, 0, 0, 0, 0, 0.1, 0.1, 0.1, 0.1, 0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 0.3, 0.3, 0.3, 0.3, 0.3, 0.4, 0.4, 0.4, 0.4, 0.4, 0.5, 0.5, 0.5, 0.5, 0.5, 0, 0, 0, 0, 0, 0.1, 0.1, 0.1, 0.1, 0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 0.3, 0.3, 0.3, 0.3, 0.3, 0.4, 0.4, 0.4, 0.4, 0.4, 0.5, 0.5, 0.5, 0.5, 0.5, 0, 1, 2, 2.5, 0, 1, 2, 2.5, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 4, 0, 4, 0, 4, 0, 4), Delta = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5), sigma = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4, 4, 7, 7, 7, 7, 7, 7, 7, 7, 4, 4, 4, 4, 8, 8, 8, 8, 10, 10, 10, 10, 14, 14, 14, 14, 35, 35, 30, 50, 6, 6, 9, 9), n1 = c(10, 15, 20, 25, 30, 10, 15, 20, 25, 30, 10, 15, 20, 25, 30, 10, 15, 20, 25, 30, 10, 15, 20, 25, 30, 10, 15, 20, 25, 30, 10, 15, 20, 25, 30, 10, 15, 20, 25, 30, 10, 15, 20, 25, 30, 10, 15, 20, 25, 30, 10, 15, 20, 25, 30, 10, 15, 20, 25, 30, 50, 50, 50, 50, 10, 10, 10, 10, 100, 100, 100, 100, 100, 100, 100, 100, 185, 185, 185, 185, 185, 185, 185, 185, 250, 250, 250, 250, 500, 500, 500, 500, 600, 600, 600, 600, 1190, 1190, 1190, 1190), n2 = c(10, 15, 20, 25, 30, 10, 15, 20, 25, 30, 10, 15, 20, 25, 30, 10, 15, 20, 25, 30, 10, 15, 20, 25, 30, 10, 15, 20, 25, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 50, 50, 50, 50, 90, 90, 90, 90, 100, 100, 100, 100, 100, 100, 100, 100, 10, 10, 10, 10, 100, 100, 100, 100, 250, 250, 250, 250, 500, 500, 500, 500, 600, 600, 600, 600, 10, 10, 10, 10), powerPASS = c(0.3909, 0.6954, 0.8558, 0.9343, 0.9709, 0.3827, 0.6784, 0.8366, 0.9178, 0.9589, 0.3591, 0.6295, 0.7807, 0.868, 0.9202, 0.3229, 0.5554, 0.6932, 0.7847, 0.8491, 0.2782, 0.4653, 0.5831, 0.6717, 0.7426, 0.2297, 0.3697, 0.4621, 0.5388, 0.606, 0.7037, 0.8576, 0.9232, 0.9545, 0.9709, 0.6865, 0.8385, 0.906, 0.94, 0.9589, 0.637, 0.7826, 0.8544, 0.895, 0.9202, 0.5619, 0.695, 0.7694, 0.8168, 0.8491, 0.4706, 0.5847, 0.6561, 0.7059, 0.7426, 0.3736, 0.4634, 0.5247, 0.5705, 0.606, 0.9624, 0.7985, 0.3433, 0.1529, 0.8179, 0.7035, 0.436, 0.2982, 0.9828, 0.9151, 0.6438, 0.2617, 0.9083, 0.7492, 0.3744, 0.0929, 0.8457, 0.7303, 0.4546, 0.19, 0.9824, 0.9142, 0.6423, 0.261, 0.9671, 0.8452, 0.4616, 0.1129, 0.9714, 0.8552, 0.4733, 0.1161, 0.1179, 0.0169, 0.7856, 0.0268, 0.2343, 0.0277, 0.0836, 0.0315), powerSAS = c(0.39094, 0.69541, 0.8558, 0.93426, 0.97092, 0.38272, 0.67836, 0.83661, 0.91781, 0.95889, 0.35908, 0.62954, 0.78069, 0.86796, 0.92017, 0.32287, 0.5554, 0.69322, 0.78471, 0.84909, 0.2782, 0.46531, 0.58306, 0.67174, 0.7426, 0.2297, 0.36967, 0.46208, 0.53883, 0.606, 0.70373, 0.85764, 0.92324, 0.95452, 0.97092, 0.68648, 0.83846, 0.90604, 0.94003, 0.95889, 0.63703, 0.78255, 0.85439, 0.89501, 0.92017, 0.56192, 0.69503, 0.76944, 0.81676, 0.84909, 0.47061, 0.5847, 0.6561, 0.70594, 0.7426, 0.37365, 0.46343, 0.52475, 0.57052, 0.606, 0.96239, 0.79849, 0.34329, 0.15288, 0.81791, 0.70346, 0.43595, 0.29819, 0.98278, 0.91512, 0.64376, 0.26169, 0.90831, 0.74924, 0.37443, 0.0929, 0.84568, 0.73025, 0.45459, 0.19001, 0.9824, 0.91417, 0.64226, 0.26103, 0.96713, 0.84523, 0.46161, 0.11288, 0.97112, 0.85433, 0.47184, 0.11536, 0.11547, 0.01658, 0.78512, 0.02642, 0.23194, 0.02739, 0.08256, 0.03115 )), .Names = c("alpha", "delta0", "Delta", "sigma", "n1", "n2", "powerPASS", "powerSAS"), class = "data.frame", row.names = c(NA, -100L)) SAS <- SAS[,-7] ## ----echo=FALSE--------------------------------------------------------------- knitr::kable(SAS, row.names = 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] ) } ## ----------------------------------------------------------------------------- identical(round(power,5), SAS$powerSAS) ## ----------------------------------------------------------------------------- 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) ## ----------------------------------------------------------------------------- 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) } ## ----echo=FALSE, fig.width=8, fig.height=4------------------------------------ oldpar <- par(mar=c(4,4,0.4,0.4)) layout(t(c(1,2))) sigma <- seq(65,69,len=100) n1 <- 1000; n2 <- 1000 plot(sigma, powerTOST(0.05, 0, 5, sigma, n1, n2), type="l", lwd=2, xlab=expression(sigma), ylab="power") y <- sapply(sigma, function(sigma) ipowerTOST(0.05, 0, 5, sigma, n1, n2)) lines(sigma, y, col="blue", lwd=2) y <- sapply(sigma, function(sigma) powerTOST(0.05, 0, 5, sigma, n1, n2, algo=1)) lines(sigma, y, col="red", lwd=2) abline(h=0, col="green", lty=2) legend("topright", c("powen4 - 1", "powen4 - 2", "ipowen4"), lty=c(1,1,1), col=c("red", "black", "blue")) sigma <- seq(15,69,len=100) plot(sigma, powerTOST(0.05, 0, 5, sigma, n1, n2), type="l", lwd=2, xlab=expression(sigma), ylab="power") y <- sapply(sigma, function(sigma) ipowerTOST(0.05, 0, 5, sigma, n1, n2)) lines(sigma, y, col="blue", lwd=2) y <- sapply(sigma, function(sigma) powerTOST(0.05, 0, 5, sigma, n1, n2, algo=1)) lines(sigma, y, col="red", lwd=2) abline(h=0, col="green", lty=2) legend("topright", c("powen4 - 1", "powen4 - 2", "ipowen4"), lty=c(1,1,1), col=c("red", "black", "blue")) par(oldpar) ## ----echo=FALSE, fig.width=4, fig.height=4------------------------------------ oldpar <- par(mar=c(4,4,0.4,0.4)) n1 <- n2 <- 720 sigma <- seq(56,57,len=100) plot(sigma, powerTOST(0.05, 0, 5, sigma, n1, n2), type="l", lwd=2, xlab=expression(sigma), ylab="power") y <- sapply(sigma, function(sigma) ipowerTOST(0.05, 0, 5, sigma, n1, n2)) lines(sigma, y, col="blue", lwd=2) y <- sapply(sigma, function(sigma) powerTOST(0.05, 0, 5, sigma, n1, n2, algo=1)) lines(sigma, y, col="red", lwd=2) legend("topright", c("powen4 - 1", "powen4 - 2", "ipowen4"), lty=c(1,1,1), col=c("red", "black", "blue")) par(oldpar) ## ----echo=FALSE, fig.width=8, fig.height=4------------------------------------ oldpar <- par(mar=c(4, 4, 0.2, 0.2)) layout(t(c(1,2))) n1 <- n2 <- 700 sigma <- seq(35,45,len=100) plot(sigma, powerTOST(0.01, 1, 5, sigma, n1, n2), type="l", lwd=2, xlab=expression(sigma), ylab="power") y <- sapply(sigma, function(sigma) ipowerTOST(0.01, 1, 5, sigma, n1, n2)) lines(sigma, y, col="blue", lwd=2) y <- sapply(sigma, function(sigma) powerTOST(0.05, 0, 5, sigma, n1, n2, algo=1)) lines(sigma, y, col="red", lwd=2) legend("topright", c("powen4 - 1", "powen4 - 2", "ipowen4"), lty=c(1,1,1), col=c("red", "black", "blue")) n1 <- n2 <- 700 sigma <- seq(38.5,39,len=100) plot(sigma, powerTOST(0.01, 1, 5, sigma, n1, n2), type="l", lwd=2, xlab=expression(sigma), ylab="power") y <- sapply(sigma, function(sigma) ipowerTOST(0.01, 1, 5, sigma, n1, n2)) lines(sigma, y, col="blue", lwd=2) y <- sapply(sigma, function(sigma) powerTOST(0.05, 0, 5, sigma, n1, n2, algo=1)) lines(sigma, y, col="red", lwd=2) legend("topright", c("powen4 - 1", "powen4 - 2", "ipowen4"), lty=c(1,1,1), col=c("red", "black", "blue")) par(oldpar) ## ----echo=FALSE, fig.width=4, fig.height=4------------------------------------ oldpar <- par(mar=c(4, 4, 0.7, 0.2)) n1 <- n2 <- 600 sigma <- seq(30,36,len=100) plot(sigma, powerTOST(0.005, 0, 5, sigma, n1, n2), type="l", lwd=2, xlab=expression(sigma), ylab="power") y <- sapply(sigma, function(sigma) ipowerTOST(0.005, 0, 5, sigma, n1, n2)) lines(sigma, y, pch=19, col="blue", lwd=2) y <- sapply(sigma, function(sigma) powerTOST(0.05, 0, 5, sigma, n1, n2, algo=1)) lines(sigma, y, col="red", lwd=2) legend("topright", c("powen4 - 1", "powen4 - 2", "ipowen4"), lty=c(1,1,1), col=c("red", "black", "blue")) par(oldpar) ## ----------------------------------------------------------------------------- alpha <- 0.05; delta0 <- 0; Delta <- 5 sigma <- 110 n1 <- n2 <- 2500 powerTOST(alpha, delta0, Delta, sigma, n1, n2) ipowerTOST(alpha, delta0, Delta, sigma, n1, n2) ## ----------------------------------------------------------------------------- sigma <- 152 n1 <- n2 <- 5000 powerTOST(alpha, delta0, Delta, sigma, n1, n2) ipowerTOST(alpha, delta0, Delta, sigma, n1, n2) ## ----------------------------------------------------------------------------- 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)) ## ----------------------------------------------------------------------------- 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)) ## ----------------------------------------------------------------------------- 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) ## ----------------------------------------------------------------------------- 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) ## ----------------------------------------------------------------------------- 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) ## ----------------------------------------------------------------------------- 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)