Sampling distribution of eta^2, omega^2-hat, and epsilon^2

sample.eta2.R

## 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)

utility_functions.R

## 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)
  }
}