title "共分散構造分析(csm,SEM)の必要サンプル数".
* 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=100000.
data list free/ alpha rmsea0 rmseaa powd df.
begin data
0.05 0.08 0.05 .80 4
0.05 0.05 0.01 .80 4
0.05 0.08 0.05 .80 6
0.05 0.05 0.01 .80 6
end data.
*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.
*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>=.001).
  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(f8.0) alpha rmsea0 rmseaa powa powd (f7.3).
report format=list/var= rmsea0 'null' 'hyp' 'RMSEA' rmseaa  'alt' 'hyp' 'RMSEA' powd 'desired' 'power' powa 'aquired' 'power' alpha df minn '必要' 'サンプル数'.
execute.