Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Sampling distribution of eta^2, omega^2-hat, and epsilon^2
## Explore the sampling distribution of eta^2, omega^2-hat, and epsilon^2
## Richard D. Morey
## Sept 30, 2017
## Settings
## Number of participants in group
N = 10
## Number of groups
J = 3
## True effect size (prop. of var.)
omega2 = .2
#### No need to change anything below
## Define utility functions (from github gist)
source('https://gist.githubusercontent.com/richarddmorey/ea32db6a6c2419eaff0e1a5186e986c0/raw/1a78da1e938b96be400f5495156a48ba3267b378/utility_functions.R')
## Compute other needed/useful quantities
df1 = J - 1
df2 = J * (N - 1)
d.oh2 = (df2 + 1)/df1
d.eps = df2 / df1
ncp = omega2 / (1 - omega2) * N * J
exp.eta2 = compute.moment(deta2, df1, df2, ncp, min = 0)
exp2.eta2 = compute.moment(deta2, df1, df2, ncp, moment = 2, min = 0)
var.eta2 = exp2.eta2 - exp.eta2^2
sd.eta2 = sqrt(var.eta2)
exp.oh2 = compute.moment(doh2, df1, df2, ncp, min = -1/d.oh2)
exp2.oh2 = compute.moment(doh2, df1, df2, ncp, moment = 2, min = -1/d.oh2)
var.oh2 = exp2.oh2 - exp.oh2^2
sd.oh2 = sqrt(var.oh2)
exp.eps2 = compute.moment(deps2, df1, df2, ncp, min = -1/d.eps)
exp2.eps2 = compute.moment(deps2, df1, df2, ncp, moment = 2, min = -1/d.eps)
var.eps2 = exp2.eps2 - exp.eps2^2
sd.eps2 = sqrt(var.eps2)
## Output ###########
### eta^2
## Expected value
exp.eta2
## Bias when estimating omega2
exp.eta2 - omega2
## % bias
(exp.eta2 - omega2) / omega2 * 100
## Standard deviation
sd.eta2
## RMSE
sqrt((var.eta2 + (exp.eta2 - omega2)^2))
## Selected quantiles
qeta2(c(.025, .5, .975), df1, df2, ncp)
## Probability that eta^2 > omega^2
## (i.e., over-estimation probability)
1 - peta2(omega2, df1, df2, ncp)
### omega^2-hat
## Expected value
exp.oh2
## Bias when estimating omega2
exp.oh2 - omega2
## % bias
(exp.oh2 - omega2) / omega2 * 100
## Standard deviation
sd.oh2
## RMSE
sqrt((var.oh2 + (exp.oh2 - omega2)^2))
## Selected quantiles
qoh2(c(.025, .5, .975), df1, df2, ncp)
## Probability that hat(omega^2) > omega^2
## (i.e., over-estimation probability)
1 - poh2(omega2, df1, df2, ncp)
## Probability that hat(omega^2) < 0
## (i.e., meaningless estimate)
poh2(0, df1, df2, ncp)
### epsilon^2
## Expected value
exp.eps2
## Bias when estimating omega2
exp.eps2 - omega2
## % bias
(exp.eps2 - omega2) / omega2 * 100
## Standard deviation
sd.eps2
## RMSE
sqrt((var.eps2 + (exp.eps2 - omega2)^2))
## Selected quantiles
qeps2(c(.025, .5, .975), df1, df2, ncp)
## Probability that eps^2 > omega^2
## (i.e., over-estimation probability)
1 - peps2(omega2, df1, df2, ncp)
## Probability that hat(omega^2) < 0
## (i.e., meaningless estimate)
peps2(0, df1, df2, ncp)
## Plot the sampling distributions
## blue = eta^2
## red = omega^2-hat
## green, dashed = epsilon^2
curve(deta2(x, df1, df2, ncp), -1/d.eps, 1, ylab = "Density",
xlab = expression(eta^2), yaxt = "n", lwd = 2,
col = "blue", n = 501)
curve(doh2(x, df1, df2, ncp), -1/d.oh2, 1,
lwd = 2, add = TRUE, col = "red", n = 501)
curve(deps2(x, df1, df2, ncp), -1/d.eps, 1,
lwd = 2, add = TRUE, col = "darkgreen",
lty=2, n = 501)
## Check with a simulation
hist(reta2(10000, df1, df2, ncp), freq=FALSE,
col = rgb(0,0,1,.3), border = NA, add = TRUE)
## Check with a simulation
hist(roh2(10000, df1, df2, ncp), freq=FALSE,
col = rgb(1,0,0,.3), border = NA, add = TRUE)
## Check with a simulation
hist(reps2(10000, df1, df2, ncp), freq=FALSE,
col = rgb(0,1,0,.3), border = NA, add = TRUE)
## Define useful functions
compute.moment = function(FUN, df1, df2, ncp, moment = 1, min){
integrate(function(x){
x ^ moment * FUN(x, df1, df2, ncp)
}, min, 1)[[1]]
}
## Density function of eta^2
deta2 = function(e2, df1, df2, ncp){
d = df2 / df1
f = e2 / (1 - e2) * d
df(f, df1, df2, ncp) * d * ( 1 / (1 - e2) + e2 / (1 - e2)^2 )
}
## CDF of eta^2
peta2 = function(e2, df1, df2, ncp){
d = df2 / df1
f = e2 / (1 - e2) * d
pf(f, df1, df2, ncp)
}
## Quantile function of eta^2
qeta2 = function(p, df1, df2, ncp){
d = df2 / df1
f = qf(p, df1, df2, ncp)
f / (f + d)
}
## Random generation of eta^2
reta2 = function(n, df1, df2, ncp){
d = df2 / df1
f = rf(n, df1, df2, ncp)
f / (f + d)
}
## Density function of omega-hat^2
doh2 = function(oh2, df1, df2, ncp){
d = (df2 + 1)/df1
f = ( oh2 * d + 1 ) / (1 - oh2)
df(f, df1, df2, ncp) * abs( d / (1 - oh2) + (oh2 * d + 1) / (1 - oh2)^2 )
}
## CDF of omega-hat^2
poh2 = function(oh2, df1, df2, ncp){
d = (df2 + 1)/df1
f = (oh2 * d + 1) / (1 - oh2)
pf(f, df1, df2, ncp)
}
## Quantile function of omega-hat^2
qoh2 = function(p, df1, df2, ncp){
d = (df2 + 1)/df1
f = qf(p, df1, df2, ncp)
(f - 1) / ( f + d )
}
## Random generation of omega-hat^2
roh2 = function(n, df1, df2, ncp){
d = (df2 + 1)/df1
f = rf(n, df1, df2, ncp)
(f - 1) / ( f + d )
}
## Density function of eps^2
deps2 = function(eps2, df1, df2, ncp){
d = df2/df1
f = ( eps2 * d + 1 ) / (1 - eps2)
df(f, df1, df2, ncp) * abs( d / (1 - eps2) + (eps2 * d + 1) / (1 - eps2)^2 )
}
## CDF of eps^2
peps2 = function(eps2, df1, df2, ncp){
d = df2/df1
f = (eps2 * d + 1) / (1 - eps2)
pf(f, df1, df2, ncp)
}
## Quantile function of eps^2
qeps2 = function(p, df1, df2, ncp){
d = df2/df1
f = qf(p, df1, df2, ncp)
(f - 1) / ( f + d )
}
## Random generation of eps^2
reps2 = function(n, df1, df2, ncp){
d = df2/df1
f = rf(n, df1, df2, ncp)
(f - 1) / ( f + d )
}
### Confidence interval
steigerCI2.omega2 = function(fstat, df1, df2, conf.level = .95){
alpha = 1 - conf.level
c(steigerCI2.omega2.one(fstat, df1, df2, 1 - alpha/2),
steigerCI2.omega2.one(fstat, df1, df2, alpha/2))
}
steigerCI2.omega2.one<-function(fstat, df1, df2, alpha){
if(fstat < qf(alpha, df1, df2)){
return(0)
}else{
return(optimize(function(o2){
(fstat - qf(alpha, df1, df2, ncp = (df1 + df2 + 1) * o2 / (1 - o2)))^2
}, interval = c(0,1))$minimum)
}
}
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.