R/OwensQ.R

In PowerTOST: Power and Sample Size for (Bio)Equivalence Studies

#------------------------------------------------------------------------------
# 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