*chisize.sps ------------------------------------.
* "χ2検定 必要サンプル数（表 指定型マクロ）".
*  effect size index .
* small w=.10 : medium w=.30: large w=.50  Cohen(1988).
* Cohen, J. (1988). Statistical power analysis for the behavioral sciences (2nd edition).
*          Hillsdale, NJ:Erlbaum.
define chisize (rk=!charend('/')
              /value=!charend('/')
              /row=!default(0)!charend('/')
              /n=!default(0)!charend('/')
              /alpha=!default(.05)!charend('/') 
              /power=!default(.80)!charend('/')
              ).
preserve.
set printback=none.
set mxloops=10000000.
* row=-1 proportion is deterrmined by value.
* n =-1 ntotal is sum of row or value.

!let !ncell=x .  
!do !i !in (!value) .
  !let !ncell=!concat(!ncell,x).
!doend.

!let !ncell=!length(!substr(!ncell,2)).
!let !dend=!concat(d,!ncell).
!let !r=!head(!rk).
!let !k=!tail(!rk).
!let !trowend=!concat(trow,!r).
!let !rowend=!concat(row,!r).
!let !form=!concat(!r,f5.3).

!let !nrow=x .  
!let !first=99 .

!if (!row=!null) !then.
  !let !row0='0'.
  !else.
  !let !row0=!row.
!ifend.
!do !i !in (!row0) .
  !let !nrow=!concat(!nrow,x).
  !if (!length(!nrow)=2) !then.
       !let !first=!i.
  !ifend. 
!doend.
!let !need=9.

exec.

!if (!first='-') !then .
   !let !first='-1' .
   !let !need=1.
!ifend.

!if (!length(!nrow)>1.5) !then.
  !let !nrow=!length(!substr(!nrow,2)).
!ifend.

!if (!nrow=1) !then.
   !let !need=1.
   !else .
   !let !row0=!row .
!ifend.

!if (!need=1) !then.
  !let !row0=!concat('row1 to ',!rowend).
!ifend.

*.
input program.
    compute vdummy=1.
   end case.
end file.
end input program.
*.

compute #krow=1.
if (!nrow=1) #krow=0.
if (!nrow=2) #krow=!first.
compute powd=!power.
compute alpha=!alpha.
compute r=!r.
compute k=!k.
compute ntotal=0.
*.
vector d(!ncell).
do repeat x=d1 to !dend /y=!value.
  compute x=y.
end repeat.
vector row(!r).
vector trow(!r).
vector col(!k).
loop #j=1 to k.
 compute col(#j)=0.
end loop.
loop #i=1 to r.
 compute row(#i)=0.
 compute trow(#i)=0.
end loop.
*.
compute #kk=0.
loop #i=1 to r.
   loop #j=1 to k.
     compute #kk=#kk+1.
     compute row(#i)=row(#i)+d(#kk).
   end loop.
end loop.

compute ntotal=0.
!if (!head(!n)=!null) !then.
   !let !nn=0.
   !else.
   !let !nn=!n.
!ifend.
if (!nn>0) ntotal=!nn.
*.

do if (!need=9).
  do repeat x=trow1 to !trowend /z=!unquote(!row0).
    compute x=z.
  end repeat.
end if.

do if (#krow=-1).
  do repeat x=trow1 to !trowend /y=row1 to !rowend .
    compute x=y.
  end repeat.
end if.

do if (#krow=0).
  do repeat x=trow1 to !trowend .
    compute x=1.
  end repeat.
end if.

do if (!nn=-1).
   compute ntotal=sum(trow1 to !trowend).
end if.
compute #rsum=sum(trow1 to !trowend).


*.
compute #kk=0.
loop #i=1 to r.
   loop #j=1 to k.
     compute #kk=#kk+1.
     compute d(#kk)=d(#kk)/row(#i)*trow(#i)/#rsum.
      compute col(#j)=col(#j)+d(#kk).
   end loop.
   compute row(#i)=trow(#i)/#rsum.
end loop.
 
compute #ef=0.
compute #kk=0.
loop #i=1 to r.
   loop #j=1 to k.
     compute #kk=#kk+1.
     compute #rowcol=row(#i)*col(#j).
     compute #ef=#ef+(d(#kk)-#rowcol)**2/#rowcol.
   end loop.
   compute row(#i)=trow(#i)/#rsum.
end loop.
compute esize=sqrt(#ef).
compute w2=#ef.
compute dfh=(r-1)*(k-1).
do if (ntotal >0).
    compute lambda=w2*ntotal.
    compute c_alpha=idf.chisq(1-alpha,dfh).
    compute powa=1-ncdf.chisq(c_alpha,dfh,lambda).
    compute xalpha=1-cdf.chisq(lambda,dfh).
*    write / /'                 power   alpha   dfh  lambda  ntotal   w'
          /'sample power 'powa alpha dfh lambda ntotal esize
          /'sample alpha ' xalpha .
end if.

!if (!nn<>0) !then.
format powa alpha xalpha esize (f8.3) dfh ntotal (f4.0).
report format=list/var=
 powa 'observed' 'power'
 alpha 'desired' 'alpha'
 dfh 'df' (4)
 ntotal 'sample' 'size' (8)
 xalpha 'sample' 'alpha' (6)
 lambda 'lambda' (8)
 esize  'w' 'effect' 'size' (6)
 /title '想定したサンプル数に基づく統計量'
.
!ifend.

*.
write /'proportin (row)'.
write /row1 to !rowend ( !form //).

*initialize values.
compute powa=0.0.
compute n=0.
compute w2=esize**2.
*compute dfh=df.
loop  if  (powa<=powd).
  compute n=n +100.
  compute lambda=w2*n.
  compute c_alpha=idf.chisq(1-alpha,dfh).
  compute powa=1-ncdf.chisq(c_alpha,dfh,lambda).
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 lambda=w2*newn.
*compute power .
  compute c_alpha=idf.chisq(1-alpha,dfh).
  compute powa=1-ncdf.chisq(c_alpha,dfh,lambda).
  compute powdiff=abs(powa-powd).
  do if (powa<powd).
    compute  dir=1.
    else .
    compute dir= -1.
  end if.
end loop.
compute minn=newn.
format esize powd powa alpha (f5.3) dfh(f8.0).

report format=list/var= 
     esize 'effect' 'size ' '(w)' 
     alpha 'alpha' 'level'
     powd 'desired' 'power'  
     powa 'acquired' 'power' 
     dfh 'degree' 'of' 'freedom'
     minn 'minimum' 'sample' 'size' 'total'
     /title 'χ2検定　必要サンプル数'
 .
restore.
!enddefine.

chisize rk=2 2/value=55 45 90 110 /row=150 400/n=-1.
exec.