title "因子分析のの必要サンプル数 （cmsの応用）".
* MacCallum,R.C.,Browne,M.W.,and Sugawara,H.M.(1996) Power analysis and
* determination of sample size for covariance structure modeling.
* Psychologica Mesthods, 1(2) 130-149.3..
*http://quantrm2.psy.ohio-state.edu:80/MacCallum/power.htm にプログラムがある.
set mxloops=10000.
data list free/ alpha rmsea0 rmseaa powd vn.
begin data
0.05 0.08 0.05   .80 13
0.05 0.05 0.01   .80 13
0.05 0.08 0.05   .80 20
0.05 0.05 0.01   .80 20
0.05 0.08 0.05   .80 100
0.05 0.05 0.01   .80 100
end data.
*最大因子数 自由度>0 となるもの.
compute mxf=trunc(((2*vn+1)-sqrt(8*vn+1))/2-0.000001).  
compute mxf=min(mxf,9). 
*close fit.
*rmsea0=.05 ; *null hyp rmsea ;
*rmseaa=.08 ; *alt hyp rmsea ;
*not close fit.
*rmsea0=.05 ; *null hyp rmsea ;
*rmseaa=.01 ; *alt hyp rmsea ;
*
*df=20;       *degrees of freedom.
*alpha=.05;  *alpha level.
*powd=.80;   *desired power.df.
compute m=1.

*write /'rmsea0  rmseaa    powd   powa   alpha    df      vn      m    minn' .
write /'rmsea0  rmseaa    powd   powa   alpha    df   変数数 因子数 必要サンプル数' .


loop if (m<mxf+1).
compute df=((vn-m)**2-(vn+m))/2.
*initialize values.
compute powa=0.0.
compute n=0.
do if (rmsea0<rmseaa).
  compute beta=1-alpha.
else.
  compute beta=alpha.
end if.
compute dfcinv=df.

*begin loop for finding initial level of n.
loop  if  (powa<=powd).
  compute n=n +100.
  compute ncp0=(n-1)*df*rmsea0**2.
  compute ncpa=(n-1)*df*rmseaa**2.
  compute start=df+2*ncp0.
*compute power.
*compute cinv.
   compute lambdaci=ncp0.
   compute stepcinv=start/20.
   compute prec=.00000001.
   loop icinv=0 to 1000.
     compute cval=start-stepcinv*icinv.
   end loop if ((beta)>ncdf.chisq(cval,dfcinv,lambdaci)).
   do if (icinv=0).
     loop icinv=0 to 1000.
        compute cval=start+stepcinv*icinv.
     end loop if ((beta)<ncdf.chisq(cval,dfcinv,lambdaci)).
     compute icinv=icinv*(-1).
   end if.
   compute betaa=ncdf.chisq(cval,dfcinv,lambdaci).

   *begin loop for interval halving.
   compute dircinv=icinv/abs(icinv).
   compute newncinv=cval.
   compute intvcinv=stepcinv*2.
   compute betadiff=abs(betaa-beta).
   loop if (betadiff>=prec).
     compute  intvcinv=intvcinv*.5.
     compute  newncinv=newncinv +dircinv*intvcinv*.5.
   *compute new betaa.
     compute betaa=ncdf.chisq(newncinv,dfcinv,lambdaci).
     compute betadiff=abs(betaa-beta).
     do if (betaa<beta).
       compute  dircinv=1. 
       else .
       compute dircinv= -1.
     end if.
   end loop .
   *end of cinv.
   compute cval=newncinv.
   do if (rmsea0<rmseaa) .
   compute  powa=1-ncdf.chisq(cval,df,ncpa).
   else.
   compute  powa=ncdf.chisq(cval,df,ncpa).
   end if.
end loop.
******************************************************.
*begin loop for interval halving.
compute dir=-1.
compute newn=n.
compute intv=200.
compute powdiff=powa-powd.
loop if (powdiff>=.0001).
  compute  intv=intv*.5.
  compute  newn=newn +dir*intv*.5.
*compute new power.
  compute  ncp0=(newn-1)*df*rmsea0**2.
  compute  ncpa=(newn-1)*df*rmseaa**2.
  compute start=df+2*ncp0.
*compute power .
*compute cinv ****************************.
   compute lambdaci=ncp0.
   compute stepcinv=start/20.
   compute prec=.00000001.
   loop icinv=0 to 1000.
     compute cval=start-stepcinv*icinv.
   end loop if ((beta)>ncdf.chisq(cval,dfcinv,lambdaci)).
   do if (icinv=0).
     loop icinv=0 to 1000.
        compute cval=start+stepcinv*icinv.
     end loop if ((beta)<ncdf.chisq(cval,dfcinv,lambdaci)).
     compute icinv=icinv*(-1).
   end if.
   compute betaa=ncdf.chisq(cval,dfcinv,lambdaci).

   *begin loop for interval halving.
   compute dircinv=icinv/abs(icinv).
   compute newncinv=cval.
   compute intvcinv=stepcinv*2.
   compute betadiff=abs(betaa-beta).
   loop if (betadiff>=prec).
     compute  intvcinv=intvcinv*.5.
     compute  newncinv=newncinv +dircinv*intvcinv*.5.
   *compute new betaa.
     compute betaa=ncdf.chisq(newncinv,dfcinv,lambdaci).
     compute betadiff=abs(betaa-beta).
     do if (betaa<beta).
       compute  dircinv=1. 
       else .
       compute dircinv= -1.
     end if.
   end loop .
   *end of cinv *****************************.
   compute cval=newncinv.
   do if (rmsea0<rmseaa).
     compute  powa=1-ncdf.chisq(cval,df,ncpa).
     else.
     compute  powa=ncdf.chisq(cval,df,ncpa).
   end if.
  compute powdiff=abs(powa-powd).
  do if (powa<powd).
    compute  dir=1.
    else .
    compute dir= -1.
  end if.
end loop.
compute minn=newn.
format df vn m(f8.0) alpha rmsea0 rmseaa powa powd (f7.3).
write /rmsea0 rmseaa  powd powa alpha df vn m minn .
compute m=m+1.
end loop.
execute.