Skip to content

Instantly share code, notes, and snippets.

@richarddmorey richarddmorey/for_lakens.R
Last active Oct 15, 2017

Embed
What would you like to do?
For Lakens
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")
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.