Nothing
#------------------------------------------------------------------------------
# Owen's Q-function
# author: dlabes
#-------------------------------------------------------------------------------
# a, b must be a scalar numeric
# nu, t and delta also, no vectors allowed
# This is a variant with variables transformation
# and using the fact that Integral(0,b)+Integral(b,Inf) gives p(nct)
OwensQ <- function (nu, t, delta, a=0, b)
{
if (missing(b)) stop("Upper integration limit missing.")
if (a!=0) stop("Only a==0 implemented.")
if (length(nu)>1 | length(t)>1 | length(delta)>1 | length(a)>1 | length(b)>1)
stop("Input must be scalars!")
if (nu<1) stop("nu must be >=1.")
if(a==b) return(0)
# in case of alpha>=0.5 b is infinite
# also in case of se=0 and diffm != ltheta1 or !=ltheta2
# for that case see:
# A bivariate noncentral T-distibution with applications
# Youn Min Chou
# Communications in Statistics - Theory and Methods, 21:12, 3427-3462, 1992
# DOI: 10.1080/03610929208830988
#
# in case of abs(delta)>37.62 the non-central t is evaluated via a crude
# approximation which can be rather poor for small nu
# then we return OwensQOwen() but for speed reasons only for small nu
# what is small is more or less arbitrary
if (nu < 29 && abs(delta) > 37.62) {
if (is.infinite(b)) {
# for high delta we won't use the non-central t
# thus we try numerical integration of the original density from 0 to Inf,
# remapped to 0,1
i_fun <- function(y) .Q.integrand(y / (1 - y), nu, t, delta) * 1/(1 - y)^2
return(integrate(i_fun, lower=0, upper=1, subdivisions = 1000L,
rel.tol = 1.e-8, stop.on.error = TRUE)[[1]])
} else {
return(OwensQOwen(nu, t, delta, 0, b))
}
} else {
if (is.infinite(b)) {
#we use the nct according to the Chou paper for low to moderate delta
return(suppressWarnings(pt(t, df=nu, ncp=delta)))
} else {
# We calculate Owen's Q via
# pt(t, df=nu, ncp=delta) - Integral(b,Inf)
# The Integral(b,Inf) is via transformation of the variables x=b+y/(1-y)
# remapped to an integral(0,1)
i_fun <- function(y) {
.Q.integrand(b + y/(1-y), nu, t, delta) / (1-y)^2
}
Integral01 <- integrate(i_fun, lower=0, upper=1, subdivisions = 1000L,
rel.tol = 1.e-8, stop.on.error = TRUE)
# suppress the warning w.r.t. precision of nct
return(suppressWarnings(pt(t, df=nu, ncp=delta)) - Integral01[[1]])
}
}
}
#-------------------------------------------------------------------------------
# Integrand of the definit integral in Owen's Q. Used in the call of integrate()
# Not useful alone, I think ? Leading . hides this function
# function must give a vectorized answer in respect to x
.Q.integrand <- function(x, nu, t, delta)
{ #version without for - loop, it works without
lnQconst <- -((nu/2.0)-1.0)*log(2.0) - lgamma(nu/2.0)
# what if x<0? Should here not possible, but ...
# simple x^(nu-1) doesnt work for high nu because = Inf
# and then exp( -0.5*x^2 + lnQconst )*x^(nu-1) -> NaN
# (nu-1)*log(abs(x)) is NaN if nu=1, x=0! 0*(-Inf) -> NaN
dens <- x # assures that dens=0 if x=0
x1 <- x[x!=0]
dens[x!=0] <- sign(x1)^(nu-1) *
pnorm( t*x1/sqrt(nu) - delta, mean = 0, sd = 1, log.p = FALSE) *
exp( (nu-1)*log(abs(x1)) - 0.5*x1^2 + lnQconst )
# return
dens
}
# Test cases:
# Craig Zupke's observations:
# power.TOST(0.410,FALSE,-5.97,5.97,8.5448,1,14,"parallel",TRUE) #!old call
# power.TOST(0.410,FALSE,-5.97,5.97,8.5448,1,14,"parallel","exact")
# gave an error; high b/delta
# should give: 2.335633e-07
# Jul 2012: Helmuts observation
# n=4, CV=1E-5(=se) gives power=1 (delta1=24303.3, delta2=-38811.23, R=b=15283.88
# CV=1E-6 gives power=0 ( 243033 -388112.3 R 152838.8
# CV=0 gives power=1 Inf -Inf Inf
# tval=2.919986
# for CV=1e-6: erroneous in versions pre 0.9-9. The 2. call gave =0
# OwensQ(nu=2, t= 2.919986, delta= 243033, 0, 152838.8) ==0
# OwensQ(nu=2, t=-2.919986, delta=-388112.3, 0, 152838.8) ==1
# for CV=0 - no longer staistics
# OwensQ(nu=2, t=2.919986, delta=Inf, 0, Inf) ==0
# OwensQ(nu=2, t=-2.919986, delta=-Inf, 0,Inf) ==1
#
# Helmuts cases (ver 0.9-9) Jul 2012
# sampleN.TOST(theta0=1, CV=0.02, design="2x2", print=TRUE) # Ok
# next gave an error due to 0*-Inf in .Q.integrand()
# sampleN.TOST(theta0=1, CV=0.01, design="2x2", print=TRUE)
#
# Jiri Hofmann's case(s) (V1.1-08, Jan 2014)
# power.TOST(CV=0.39, n=7000, theta0=1.24) gave power=0.3511889, correct
# power.TOST(CV=0.385, n=7000, theta0=1.24) gave power=0, correct is 0.3568784
# power.TOST(CV=0.38, n=7000, theta0=1.24) gave power=0, correct is 0.3627554
# power.TOST(CV=0.375, n=7000, theta0=1.24) gave power=0, correct is 0.3688277
# power.TOST(CV=0.37, n=7000, theta0=1.24) gave power=0.3751032, correct
#
# power.TOST(CV=0.37, n=7000, theta0=1.0) gave power=0, correct is 1
# power.TOST(CV=0.375, n=7000, theta0=1.0) gave power=0, correct is 1
# power.TOST(CV=0.38, n=7000, theta0=1.0) gave power=0, correct is 1
# power.TOST(CV=0.385, n=7000, theta0=1.0) gave power=0, correct is 1
# power.TOST(CV=0.39, n=7000, theta0=1.0) gave power=1, correct
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.