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