#This code is designed to calculate optimal alpha levels for ANOVA, authored by Joe Mudge (joe.mudge@unb.ca). #The function used to calculate optimal alphas is: optab.anova(u=NULL,v=NULL,f2=NULL,T1T2cratio=1,HaHopratio=1) #The arguments 'u' and 'v' are the numerator and denominator degrees of freedom, respectively. #The argument 'f2' is the 'Cohen's f^2' standardized critical effect size. Cohen's f^2 = variance among group means/pooled within group variance #The argument 'T1T2cratio' is the cost ratio of Type I errors relative to Type II errors. T1T2cratio is set at 1 as a default, making Type I and Type II errors equally serious. #The argument 'HaHopratio' is the prior probability of the alternate hypothesis relative to the prior probability of the null hypothesis. HaHopratio is set at 1 as a default, to not weight alpha and beta by their prior probabilities (assuming they are unknown). #This code is partially based on code modified from the R package 'pwr'(Champely 2009) beta.glm.test<-function (u = NULL, v = NULL, f2 = NULL, sig.level = 0.05, power = NULL) { if (sum(sapply(list(u, v, f2, power, sig.level), is.null)) != 1) stop("exactly one of u, v, f2, power, and sig.level must be NULL") if (!is.null(f2) && f2 < 0) stop("f2 must be positive") if (!is.null(u) && u < 1) stop("degree of freedom u for numerator must be at least 1") if (!is.null(v) && v < 1) stop("degree of freedom v for denominator must be at least 1") if (!is.null(sig.level) && !is.numeric(sig.level) || any(0 > sig.level | sig.level > 1)) stop(sQuote("sig.level"), " must be numeric in [0, 1]") if (!is.null(power) && !is.numeric(power) || any(0 > power | power > 1)) stop(sQuote("power"), " must be numeric in [0, 1]") p.body <- quote({ lambda <- f2 * (u + v + 1) pf(qf(sig.level, u, v, lower = FALSE), u, v, lambda, lower = FALSE) }) if (is.null(power)) power <- eval(p.body) else if (is.null(u)) u <- uniroot(function(u) eval(p.body) - power, c(1 + 1e-10, 100))\$root else if (is.null(v)) v <- uniroot(function(v) eval(p.body) - power, c(1 + 1e-10, 1e+05))\$root else if (is.null(f2)) f2 <- uniroot(function(f2) eval(p.body) - power, c(1e-07, 1e+07))\$root else if (is.null(sig.level)) sig.level <- uniroot(function(sig.level) eval(p.body) - power, c(1e-10, 1 - 1e-10))\$root else stop("internal error") METHOD <- "Multiple regression power calculation" 1-power } w.average.error.glm<-function (alpha=NULL,u=NULL,v=NULL,f2=NULL,T1T2cratio=1,HaHopratio=1) ((alpha*T1T2cratio+(HaHopratio*beta.glm.test(u=u,v=v,f2=f2,sig.level=alpha))))/(HaHopratio+T1T2cratio) optimize.average.error.glm<-function (f, interval, ..., lower = min(interval), upper = max(interval), maximum = FALSE, tol = .Machine\$double.eps^0.25) { if (maximum) { val <- .Internal(fmin(function(arg) -f(arg, ...), lower, upper, tol)) list(maximum = val, objective = f(val, ...)) } else { val <- .Internal(fmin(function(arg) f(arg, ...), lower, upper, tol)) f(val, ...) } } min.average.error.glm<-function (u=NULL,v=NULL,f2=NULL,T1T2cratio=1,HaHopratio=1) optimize.average.error.glm(w.average.error.glm,c(0,1),tol=0.00000000001,u=u,v=v,f2=f2,T1T2cratio=T1T2cratio,HaHopratio=HaHopratio) optimize.alpha.glm<-function (f, interval, ..., lower = min(interval), upper = max(interval), maximum = FALSE, tol = .Machine\$double.eps^0.25) { if (maximum) { val <- .Internal(fmin(function(arg) -f(arg, ...), lower, upper, tol)) list(maximum = val, objective = f(val, ...)) } else { val <- .Internal(fmin(function(arg) f(arg, ...), lower, upper, tol)) val } } alpha.glm<-function (u=NULL,v=NULL,f2=NULL,T1T2cratio=1,HaHopratio=1) optimize.alpha.glm(w.average.error.glm,c(0,1),tol=0.000000000001,u=u,v=v,f2=f2,T1T2cratio=T1T2cratio,HaHopratio=HaHopratio) beta.glm<-function (u=NULL,v=NULL,f2=NULL,T1T2cratio=1,HaHopratio=1) ((T1T2cratio+HaHopratio)*min.average.error.glm(u=u,v=v,f2=f2,T1T2cratio=T1T2cratio,HaHopratio=HaHopratio)-T1T2cratio*alpha.glm(u=u,v=v,f2=f2,T1T2cratio=T1T2cratio,HaHopratio=HaHopratio))/HaHopratio optab.anova<-function (u=NULL,v=NULL,f2=NULL,T1T2cratio=1,HaHopratio=1) { list("output"=t(data.frame("numerator degrees of freedom"=u,"denominator degrees of freedom"=v,"Cohen's f squared effect size"=f2,"Type I/II error cost ratio"=T1T2cratio,"Ha/Ho prior probability ratio"=HaHopratio,"average of alpha and beta"=(alpha.glm(u=u,v=v,f2=f2,T1T2cratio=T1T2cratio,HaHopratio=HaHopratio)+HaHopratio*beta.glm(u=u,v=v,f2=f2,T1T2cratio=T1T2cratio,HaHopratio=HaHopratio))/(1+HaHopratio),"cost-weighted average of alpha and beta"=min.average.error.glm(u=u,v=v,f2=f2,T1T2cratio=T1T2cratio,HaHopratio=HaHopratio), "optimal alpha"=alpha.glm(u=u,v=v,f2=f2,T1T2cratio=T1T2cratio,HaHopratio=HaHopratio),"optimal beta"=beta.glm(u=u,v=v,f2=f2,T1T2cratio=T1T2cratio,HaHopratio=HaHopratio),row.names="values")) )} #This code is designed to calculate optimal alpha levels for ANOVA, authored by Joe Mudge (joe.mudge@unb.ca). #The function used to calculate optimal alphas is: optab.anova(u=NULL,v=NULL,f2=NULL,T1T2cratio=1,HaHopratio=1) #The arguments 'u' and 'v' are the numerator and denominator degrees of freedom, respectively. #The argument 'f2' is the 'Cohen's f^2' standardized critical effect size. Cohen's f^2 = variance among group means/pooled within group variance #The argument 'T1T2cratio' is the cost ratio of Type I errors relative to Type II errors. T1T2cratio is set at 1 as a default, making Type I and Type II errors equally serious. #The argument 'HaHopratio' is the prior probability of the alternate hypothesis relative to the prior probability of the null hypothesis. HaHopratio is set at 1 as a default, to not weight alpha and beta by their prior probabilities (assuming they are unknown). #This code is partially based on code modified from the R package 'pwr'(Champely 2009)