Compute p value for equivalence test

eqt_pvalue.R

data(attitude)

library(cocor)
## EQtest at alpha = .05 needs a 90% CI
cocor(~ rating + complaints | rating + learning, data = attitude, alternative = "less", 
      null.value = .2, test = "zou2007", conf.level = .9)

zou_ci_func = function(conf.level, ...)
{
  cocor::cocor(..., test = "zou2007", conf.level = conf.level)@zou2007$conf.int
}

eqt_pvalue = function(bounds, ci_func, opt_tol_scale = 1, ...){
  default_tol = eval(formals(stats::optimize)[["tol"]])
  opt_tol = default_tol * opt_tol_scale
  if(!is.numeric(bounds))
    stop("Bounds must be numeric.")
  switch(as.character(length(bounds)),
         "1" = c(-1,1)*abs(bounds),
         "2" = sort(bounds),
         stop("Length of bounds must be 1 or 2.")
  ) -> bounds
  estimate = ci_func(conf.level = 0, ...)[1]
  inside = estimate > bounds[1] & estimate < bounds[2]
  if(inside){
    opt_cl = optimize(function(conf.level,...){
      ci = ci_func(conf.level,...)
      if(ci[1] < bounds[1] | ci[2] > bounds[2]){
        return(.Machine$double.xmax)
      }else{
        min( (ci - bounds)^2 )
      }
    },interval = c(0,1), tol = opt_tol, ...)$minimum
    p_value = (1-opt_cl)/2
  }else{
    which_bound = ( estimate < bounds[1] ) + 1
    opt_cl = optimize(function(conf.level,...){
      (ci_func(conf.level,...)[ which_bound ] - bounds[ 3 - which_bound ])^2
    },interval = c(0,1), tol = opt_tol,...)$minimum
    p_value = (1+opt_cl)/2
  }
  return(p_value)
}

## Single equivalence bound gets turned to (-x, x)
eqt_pvalue(bounds = .4, ci_func = zou_ci_func,
           formula = ~ rating + complaints | rating + learning, data = attitude) 

## You can define any bounds, though
eqt_pvalue(bounds = c(.1,.4), ci_func = zou_ci_func,
           formula = ~ rating + complaints | rating + learning, data = attitude) 

## Should yield .05 when the bounds are the 90% CI
## First get CI
ci90 = zou_ci_func(.9, formula = ~ rating + complaints | rating + learning, data = attitude)
eqt_pvalue(bounds = ci90, ci_func = zou_ci_func,
           formula = ~ rating + complaints | rating + learning, data = attitude) 
## ...within .02%. You could decrease the tolerance to get closer

## Should be closer to .05
eqt_pvalue(bounds = ci90, ci_func = zou_ci_func, opt_tol_scale = .1,
           formula = ~ rating + complaints | rating + learning, data = attitude) 
## ...now within .004%. 

## If you want to ignore a bound (e.g., at 0) you can just set it
## to infinite. This is useful when you have bounded parameters.
## As an example, this should yield the one-sided p value for the main test 
## (oddly, a p value isn't returned by the function)
eqt_pvalue(bounds = c(0, Inf), ci_func = zou_ci_func, opt_tol_scale = .1,
           formula = ~ rating + complaints | rating + learning, data = attitude)