For Lakens

for_lakens.R

find.ncp = Vectorize(function(a=.05, b=.2, df=1){
  cr = qchisq(1-a,df)
  q = optimize(function(q){
    (pchisq(cr,df,ncp=q/(1 - q)) - b)^2
  },interval = c(0,1))$minimum
  q / (1 - q)
},"b")

find.eb = Vectorize(function(a = 0.05, b = .2, df = 1){
  ncp = find.ncp(a, b, df)
  qchisq(a, 1, ncp)
},"b")

find.p = Vectorize(function(a = 0.05, b = .2, df = 1){
  1 - pchisq(find.eb(a,b,df),1)
},"a")

alpha = .05
beta = .8
crit = qchisq(1-alpha,1)
ncp = find.ncp(a=alpha, b=1-beta)
eb = find.eb(alpha, 1-beta)
p0 = find.p(alpha, 1-beta)


xx = seq(0, 50, len = 500)
plot(xx, dchisq(xx, 1, ncp), ty='l', col="blue",lwd=2, lty=2,
     axes=FALSE,xaxs='i',yaxs='i', xpd = TRUE,
     ylab = "Density", xlab = "Samp. dist. of Z^2")
abline(v=crit, lty=2, col="red")
abline(h = 0, col = rgb(0,0,0,.4))
lines(xx, dchisq(xx, 1), col="red",lwd=2, lty=2)

N = seq(10,100)

pchisq(find.eb(alpha, 1-pow),1,find.ncp(alpha,1-pow))

pchisq(find.eb(alpha, ),1)

################
find.ncp.f = Vectorize(function(a=.05, b=.2, df1=1,df2=1){
  cr = qf(1-a,df1,df2)
  q = optimize(function(q){
    (pf(cr,df1,df2,ncp=q/(1-q)) - b)^2
  },interval = c(0,1))$minimum
  q / (1 - q)
},"df2")

find.eb.f = Vectorize(function(a = 0.05, b = .2, df1 = 1,df2 = 1){
  ncp = find.ncp.f(a, b, df1, df2)
  qf(a, df1, df2, ncp)
},"df2")

find.p.f = Vectorize(function(a = 0.05, b = .2, df1 = 1, df2 = 1){
  1 - pf(find.eb.f(a,b,df1, df2),df1,df2)
},"df2")




alpha = .05
beta = .2
N = 6:500
J = 2
df1 = J-1
df2 = J*(N-1)
crit = qf(1-alpha,df1,df2)
ncp = find.ncp.f(a=alpha, b=beta, df1, df2)
eb = find.eb.f(alpha, beta, df1, df2)
p0 = find.p.f(alpha, beta, df1, df2)

# check results
# should be alpha
1 - pf(crit, df1, df2)
# Should be power, except for optimization errors
1 - pf(crit, df1, df2, ncp)
# should be alpha, except for optimization errors
pf(eb, df1, df2, ncp)


plot(N, p0,ty='l',log="x",ylim=c(0,.3),las=1,ylab="Critical p",main="One way ANOVA")
mtext("alpha=.05, beta = .2",3,.1,adj=1)
text(max(N),max(p0),paste0("J=",J), adj=c(1,1))

J = 3
df1 = J-1
df2 = J*(N-1)
p0 = find.p.f(alpha, beta, df1, df2)
lines(N, p0, col = "red")
text(max(N),max(p0),paste0("J=",J), adj=c(1,1), col="red")

J = 4
df1 = J-1
df2 = J*(N-1)
p0 = find.p.f(alpha, beta, df1, df2)
lines(N, p0, col = "blue")
text(max(N),max(p0),paste0("J=",J), adj=c(1,1), col="blue")


J = 5
df1 = J-1
df2 = J*(N-1)
crit = qf(1-alpha,df1,df2)
ncp = find.ncp.f(a=alpha, b=beta, df1, df2)
eb = find.eb.f(alpha, beta, df1, df2)
p0 = find.p.f(alpha, beta, df1, df2)

# check results
# should be alpha
1 - pf(crit, df1, df2)
# Should be power, except for optimization errors
1 - pf(crit, df1, df2, ncp)
# should be alpha, except for optimization errors
pf(eb, df1, df2, ncp)

lines(N, p0, col = "darkgreen")
text(max(N),max(p0),paste0("J=",J), adj=c(1,1), col="darkgreen")