power_analysis_rmsea.R

power_analysis_rmsea.R

#Power analysis for CSM

alpha <- 0.05 #alpha level
d <- 29 #degrees of freedom
n <- 185 #sample size
rmsea0 <- 0.08 #null hypothesized RMSEA
rmseaa <- 0.05 #alternative hypothesized RMSEA

#Code below this point need not be changed by user
ncp0 <- (n-1)*d*rmsea0^2
ncpa <- (n-1)*d*rmseaa^2

#Compute power
if(rmsea0<rmseaa) {
    cval <- qchisq(alpha,d,ncp=ncp0,lower.tail=F)
    pow <- pchisq(cval,d,ncp=ncpa,lower.tail=F)
}
if(rmsea0>rmseaa) {
    cval <- qchisq(1-alpha,d,ncp=ncp0,lower.tail=F)
    pow <- 1-pchisq(cval,d,ncp=ncpa,lower.tail=F)
}
print(pow)

#Computation of minimum sample size for test of fit

rmsea0 <- 0.05 #null hypothesized RMSEA
rmseaa <- 0.08 #alternative hypothesized RMSEA
d <- 29 #degrees of freedom
alpha <- 0.05 #alpha level
desired <- 80 #desired power

#Code below need not be changed by user
#initialize values
pow <- 0.0
n <- 0
#begin loop for finding initial level of n
while (pow<desired) {
  n <- n+100
  ncp0 <- (n-1)*d*rmsea0^2
  ncpa <- (n-1)*d*rmseaa^2
  #compute power
  if(rmsea0<rmseaa) {
      cval <- qchisq(alpha,d,ncp=ncp0,lower.tail=F)
      pow <- pchisq(cval,d,ncp=ncpa,lower.tail=F)
  }
  else {
      cval <- qchisq(1-alpha,d,ncp=ncp0,lower.tail=F)
      pow <- 1-pchisq(cval,d,ncp=ncpa,lower.tail=F)
  }
}

#begin loop for interval halving
foo <- -1
newn <- n
interval <- 200
powdiff <- pow - desired
while (powdiff>.001) {
  interval <- interval*.5
  newn <- newn + foo*interval*.5
  ncp0 <- (newn-1)*d*rmsea0^2
  ncpa <- (newn-1)*d*rmseaa^2
  #compute power
  if(rmsea0<rmseaa) {
      cval <- qchisq(alpha,d,ncp=ncp0,lower.tail=F)
      pow <- pchisq(cval,d,ncp=ncpa,lower.tail=F)
  }
  else {
      cval <- qchisq(1-alpha,d,ncp=ncp0,lower.tail=F)
      pow <- 1-pchisq(cval,d,ncp=ncpa,lower.tail=F)
  }
  powdiff <- abs(pow-desired)
  if (pow<desired) {
    foo <- 1
  }
  if (pow>desired) {
    foo <- -1
  }
}

minn <- newn
print(minn)