Skip to content

Instantly share code, notes, and snippets.

@richarddmorey

richarddmorey/eqt_pvalue.R

Last active Jul 17, 2020
Embed
What would you like to do?
Compute p value for equivalence test
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.