SUB:数量化２類マクロ for spss
数量化２類マクロ for spss. 説明書
1996.8.9. K.Hori(hori@ec.kagawa-u.ac.jp) 
1997.01.21 改訂版 基準変数のカテゴリー数が３以上のときのバグを直す。
2000.03.17 サンプル数を 20000 まで対応するよう set mxloops=20000.
-------------------------------------------
(1)ＳＰＳＳに行列言語とマクロ言語が備わっているシステムでご利用くださ
い。
　spss for windows の行列言語は Advanced statistics にあります。マクロ
言語は Base system に含まれています。
　dos 版には何れも含まれていません。Ｍａｃ版、ＵＮＩＸ版は各自で調べて
下さい。Ｍａｃ版は古いバージョン(4)でもマクロ言語はあったようです。行
列言語があったかどうかは知りません。
---------------------------------------------
(2)出力
(a)カテゴリ数量
(b)個体数量（オプション）
(c)グループ別の個体数量の平均値、ＳＤ
(d)各変数間の（個体数量の）相関係数
(e)基準変数とアイテム間の偏相関係数
(f)相関比、重相関係数
(g)カテゴリ数量のアイテムごとの範囲（range）
(h)２類への変数投入順の固有値（相関比）の変化
(i)カレントファイルへの出力は個体数量（オプション）、各アイテムの２以
上のカテゴリに対応したダミー変数（オプション）

(i)の基準変数とダミー変数を使用してdiscriminant（windows 版ならprofess
ional のオプションが必要） を実行すれば、正準判別得点の散布図、勢力図
（値は判別分析のもの）、判別結果の要約表（ｈｉｔ表）を得ることができま
す。
---------------------------------------------
(3)指定の仕方
hayasi2 ivar=A1 A2 A3 A5/cvar=A4.


(a)ivar=説明アイテム変数リスト（ＳＰＳＳ流の省略形も使用可）
　　　　変数リストの順に投入していきます。
(b)cvar=外的基準変数（１つだけ）

         -----以上２つの指定だけで動きます。-----------

hayasi2 ivar=A1 A2 A3 A5/cvar=A4/ascore=yes/missing=yes.

(c)maxis=軸の数（既定値５または［グループ数−１］の小さい方）特にグルー
プ数が７以上で軸を５以上を必要とするなら指定する。普通は指定しない。１
６以上の指定ではエラーが生じるかもしれません。

(d)pscore= （既定値 NO）NO 以外なら個人数量を出力ウィンドウに出力（既
定値は出力しない）
(e)ascore=（既定値 NO）NO 以外なら個人数量をカレントファイル（データ窓）
に出力する。
(f)dummy=（既定値 NO）NO 以外なら説明アイテムの２以上のカテゴリに対応
したダミー変数をカレントファイル（データ窓）に出力する。変数名はCOL1,
COL2,...の続き番号。（現在のマクロのレベルではアイテムごとの番号はつけ
ることができない）
(g)missing=（既定値 NO）NO 以外なら missing のデータをlistwiseに削除す
る。
　　（これを指定すると時間がかかるでしょうが、missingデータがある場合
は指定してください。）（既定値はmissing データのチェックをしない）

 ☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆
☆既定値は たとえばpscore= そのものも入れないときの値、また ☆
☆は pscore= のあとが空白の場合の値です。                   ☆
 ☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆☆
---------------------------------------------
(3)マクロの走らせ方
(a)マクロ部分を \spsswin\ のディレクトリに保存しておく。
  ファイル名 "hayasi2.sps"
(b)次のファイルをシンタックスウィンドウに入れる
-------------------------------------------------------
include "hayasi2.sps".
hayasi2 ivar=A1 A2 A3 A5/cvar=A4/pscore=yes/missing=NO.
-------------------------------------------------------
(c)シンタックスウィンドウを上の範囲を反転させ、ｃｔｒｌ＋ｒ（コントロー
ルキーを押しながらｒを押す）。
　全範囲の指定で良ければ：　ｃｔｒｌ＋ａ、ｃｔｒｌ＋ｒ   
(d)同じセッションでは次には　include 文は必要ない（あると誤動作）。
hayasi2 ivar=A1 A2 A3 A5/cvar=A4/pscore=yes/missing=NO.
の文だけでマクロが走る。

---------------------------------------------
(4)判別分析との連動
-------------------------------------------------------------------
hayasi2 ivar=A1 a3 a2 A5/cvar=A4/dummy=yes/missing=NO.
discriminant gropus=a4(1,2)/var=col1 to col8/statistics=all/plot=all.
-------------------------------------------------------------------
数量化でdummy 出力を指定する。判別のほうでは、ダミーで出力されるcol1
 to col8
と変数をすべて指定する。これで正準判別得点の散布図、勢力図（値は判別分
析のもの）、判別結果の要約表（ｈｉｔ表）などが出力される。
　なお、今のマクロだとcol8 の８をうまく出力できないので手でいれてくだ
さい。
dummy で出力された変数名はＭａｐをみて確認して下さい。

---------------------------------------------
(5)サンプルもついていますので、不必要な部分はカットして,hayasi2.sps を作
ってspsswin に保存して下さい。

最初の実験は
(a)hayasi2.sps と実験プログラムをそれぞれ(b)の操作で別の窓に読み込む。
(b)ファイル→開く→ＳＰＳＳシンタックス→ｃｔｒｌ＋ａ、ｃｔｒｌ＋ｒ
で走らせてみて下さい。（hayasi2.spsを先に実行、その後実験プログラムを実行）
(c)実験プログラムのhayasi2 の部分をいろいろ書き換えてみる。
(d)discriminal を走らせてみる
等の実験をしてみてください。

(6)単純なチェックしかしていません。なにか問題があれば連絡を下さい。
  連絡先は堀　啓造（hori@ec.kagawa-u.ac.jp)
(7)ファイルがいくつか作成されます。ＳＰＳＳＷｉｎのディレクトリをチェ
ックしてください。
(a)'h2_tmp.sav'  オリジナルデータ(missing のチェックを指定したとき)
大きなファイルの場合削除してください。

*この下からhayasi2.sps    ------------------.
/*hayasi2.sps */
preserve.
set printback=off.
*-------------------------------------------------------------.
* 数量化２類マクロ for spss. by k.hori(hori@ec.kagawa-u.ac.) 1997/01/21
*-------------------------------------------------------------.
define hayasi2 (ivar=!charend('/') 
               /cvar=!charend('/')
               /maxis=!default(5)!charend('/')
               /pscore=!default(NO)!charend('/')
               /ascore=!default(NO)!charend('/')
               /dummy=!default(NO)!charend('/')
               /missing=!default(NO)!charend('/')
               ).
preserve.
set printback=off.
*set length=none.
*---------------------------------------------------------------------------
* Compute the number of variables
*---------------------------------------------------------------------------.
*元のデータを保存.

!let !miss2=!upcase(!missing).
!if (!miss2='') !then .
  !let !miss2='NO'.
!ifend.
!if (!miss2 <> 'NO') !then.
save outfile='h2_tmp.sav'.
count ms_= !ivar !cvar (missing).
select if ms_ = 0.
!ifend.

!let !nvs1=''.

!let !nvar = !length(!nvs1).
set mxloops=40000.
!let !dummy2=!upcase(!dummy).
!if (!dummy2='') !then .
  !let !dummy2='NO'.
!ifend.
!let !ascore2=!upcase(!ascore).
!if (!ascore2='') !then .
  !let !ascore2='NO'.
!ifend.
*****************************************************************.
matrix.
compute colname={" ","******    ","数量化２","類","******"," ",
                         "v.0.94 "}.
print colname/format=a10.
compute nprint=8.
*design の出力を1,2,3...とするため num.
compute num={1,2,3,4,5,6,7,8,9,10,11,12,13,14,15}.
compute num=t(num).
compute num2={'1','2','3','4','5','6','7','8','9','10','11','12','13','14','15'}.
get x /file=* /var=!ivar !cvar/missing=omit/sysmis=omit/names=nx.
*変数数 #nvar.
compute #nvar=ncol(nx)-1.
compute ng=nx(1,#nvar+1).
compute g=x(:,#nvar+1).
compute nx=nx(1,1:#nvar).
!let !varname=''.
* cmaxs 各itemの最大カテゴリー数.
compute cmaxs=cmax(x).
compute nsample=nrow(g).
*カテゴリの１，０展開 結果を z に.
compute z=g.
loop i = 1 to #nvar.
compute x1=x(:,i).
compute x1={num(1:cmaxs(1,i),1);x1}.
compute d=design(x1).
compute #row1=nsample+cmaxs(i).
compute d=d((cmaxs(i)+1):#row1,:).
compute d=d(:,2:cmaxs(i)).
compute z={z,d}.
end loop.
*カテゴリ総数 ctotal.
compute ctotal=ncol(z)-1.
compute x=z(:,2:(ctotal+1)).
*ダミーコードの保存（第１カテゴリは除く）.
!if (!dummy2 <> 'NO') !then .
  save x/outfile='h2_x.sav'.
!ifend. 
release z,x1,d.
*群平均値算出.
compute gr=cmax(g).
do if gr > !maxis .
  compute maxis=!maxis.
  else .
  compute maxis=gr-1.
end if.
compute x1=g.
compute x1={num(1:gr,1);x1}.
compute d=design(x1).
compute d=d((gr+1):(nsample+gr),:).
compute n=csum(d).
compute gmean=make(gr,ctotal,0).
loop i=1 to nsample.
  compute gd=make(gr,ctotal,0).
  compute gd(g(i),:)=x(i,:).
  compute gmean=gmean+gd.
end loop.
loop l=1 to gr.
compute gmean(l,:)=gmean(l,:)/n(1,l).
end loop.
compute xbar=gmean.
*---------------------------------------------------.
* Compute the pooled within groups SSCP matrix.
* Compute the beteen gropus SSCP matrix.
*---------------------------------------------------.
compute total=csum(x).
compute w=0.
loop i=1 to gr.
   compute w=w + sscp(xbar(i,:)) * n(i).
end loop.
compute w=sscp(x)-w.
compute spt=sscp(x)-sscp(total)/nsample.
compute b=spt-w.
compute sscpw=w.
compute sscpb=b.
* 本来の between within.
*compute b=b/(gr-1).
*compute w=w/(nsample-gr).
*--------------------------------.
* eigenvalue and eigen vector.
*--------------------------------.
compute k=0.
compute label={'total',num2}.
compute d=0.
compute lambda=make(#nvar,1,0).
loop i=1 to #nvar.
   compute k=k+cmaxs(1,i)-1.
   compute w=sscpw(1:k,1:k).
   compute b=sscpb(1:k,1:k).
   compute inv_w=inv(w).
   compute cw=chol(w).
   compute inv_cw=inv(cw).
   compute xx=t(inv_cw)*b*inv_cw.
   call svd(xx,evc,evl,v).
   compute evl=diag(evl).
   compute d=msum(evl/(evl+1)).
   compute lambda(i,1)={d}.
 end loop.
   print lambda/title="変数追加順と固有値（相関比）"
     /rnames=nx/cnames=label/format=f8.6.
   compute evc=inv_cw*evc.
*------------------------------------------.
*暫定のサンプル数量（判別得点）の標準偏差.
*------------------------------------------.
compute dscore=x*evc(:,1:maxis).
compute sd=sqrt((cssq(dscore)-((csum(dscore)&**2) &/nsample))/nsample).
*各カテゴリの反応数.
compute ncat=csum(x).
release x.
*固有ベクトルの調整＝＞数量化の係数.
*原点調整のための計算.
compute sumc=make(maxis,#nvar,0).
loop l=1 to maxis.
compute k=0.
loop i=1 to #nvar.
  loop j=2 to cmaxs(i).
     compute  k=k+1.
     compute sumc(l,i)=sumc(l,i)+ncat(1,k)*evc(k,l).
  end loop.
end loop.
end loop.
compute sumc=sumc &/ nsample.
compute dcoef=make(maxis,(k+#nvar),0).
*係数の原点調整と規準化.
loop l=1 to maxis.
compute k=0.  
compute k2=0.
loop i=1 to #nvar.
      compute k=k+1.
      compute dcoef(l,k)=-sumc(l,i)&/sd(1,l).
  loop j=2 to cmaxs(i).
     compute  k2=k2+1.
     compute  k=k+1.
     compute dcoef(l,k)=(-sumc(l,i)+evc(k2,l))&/sd(1,l).
  end loop.
end loop.
end loop.
compute coef=t(dcoef).
*============.
* 個体数量.
*============.
get z /file=* /var=!ivar !cvar/missing=omit/sysmis=omit.
*--------------------------------.
* x すべてのカテゴリの１、０展開.
*--------------------------------.
compute x=make(nsample,1,0).
compute #row1=0.
compute d=0.
loop i = 1 to #nvar.
   compute x1=z(:,i).
   compute maxc=cmax(x1).
   compute x1={num(1:maxc,1);x1}.
   compute d=design(x1).
   compute #row1=nsample+maxc.
   compute d=d((maxc+1):#row1,:).
   compute x={x,d}.
end loop.
release z,x1.
compute mcat=ncol(x)-1.
*.
compute x=x(:,2:(mcat+1)).
compute ncat=csum(x).
*.
compute catlab={nx(1,1),num2(1,2:cmaxs(1,1))}.
do if #nvar > 1.
loop i=2 to #nvar.
   compute catlab={catlab,nx(1,i),num2(1,2:cmaxs(1,i))}.
end loop.
end if.
compute coef={t(ncat),coef}.
compute clab={"freq",num2}.
print coef/format=f9.3/title="　　　　カテゴリ数量"
      /rnames=catlab/cnames=clab.
compute dscore=x*t(dcoef).
compute label={"group"}.
compute label={label,num2(1,1:maxis)}.
!let !pscore2=!upcase(!pscore).
!if (!pscore2='') !then.
  !let !pscore2='NO'.
!ifend.
!if (!pscore2 <> 'NO') !then .
compute dscore={g,dscore}.
print dscore/title="個体数量"/cname=label/format=f8.4.
compute dscore=dscore(:,2:(maxis+1)).
!ifend.
*================.
*  各変数の得点.
*================.
*compute l1=1.
loop l1=1 to maxis.
compute colname={" ","******","第",num2(1,l1),"相関比解","  ******"," ",ng}.
print colname/format=a8.
  compute xv=0.
  compute ek=0.
  compute vd=0.
  loop i=1 to #nvar .
    compute sk=ek+1.
    compute ek=ek+cmaxs(i).
    compute vd=block(vd,dcoef(l1,sk:ek)).
  end loop.
  compute vd=vd(2:nrow(vd),2:ncol(vd)).
  compute vd=t(vd).
  compute xv=x*vd.

*-------------------------.
*基準変数グループ別平均
*-------------------------.
compute gmean=make(1,gr,0).
compute gsd=make(1,gr,0).
loop i=1 to nsample.
  compute gd=make(1,gr,0).
  compute gd(1,g(i))=dscore(i,l1).
  compute gmean=gmean+gd.
  compute gd(1,g(i))=dscore(i,l1)*dscore(i,l1).
*  print dscore(i,l1).
*  print gd.
  compute gsd=gsd+gd.
end loop.
compute gsd=sqrt((gsd-(gmean&**2)&/(n))&/(n)).
compute gmean=gmean&/n.
compute gmsd={t(gmean),t(gsd)}.
compute clabel={"平均値","標準偏差"}.
print gmsd/format=f8.3/title="基準変数群別の平均値"/ cname=clabel
  /rname=num2.
compute gm=make(nsample,1,0).
loop i=1 to nsample.
  compute gm(i,1)=gmean(1,g(i)).
end loop.
*------------.
* 相関行列.
*------------.
  *基準変数のグループ平均値の追加.
  compute xv={xv,gm}.
  compute sumc=csum(xv).
  compute sxv=sscp(xv)-(t(sumc)*sumc)&/nsample.
  compute inv_sdxv=inv(mdiag(sqrt(diag(sxv)))).
  compute corr=inv_sdxv*sxv*inv_sdxv.
  compute nall={nx,ng}.
compute sprin=1.
compute nfac=#nvar+1.
LOOP.
compute eprin=sprin+nprint.
do if (eprin > nfac).
   compute eprin=nfac.
end if.
  compute nall2=nall(1,sprin:eprin).
  print corr(:,sprin:eprin)/
  title='相関行列'/cname=nall2/rname=nall/format=f8.3.
compute sprin=eprin+1.
END LOOP IF (nfac <= eprin).
*--------- 
* 偏相関.
*---------.
  compute inv_cor=inv(corr). 
  compute  dinv=inv(mdiag(sqrt(diag(inv_cor)))).
  compute pcorr=dinv*inv_cor*dinv.
  compute pcorr=pcorr*(-1).
  compute pcorr2=pcorr((#nvar+1),1:#nvar).
compute sprin=1.
compute nfac=#nvar.
LOOP.
compute eprin=sprin+nprint.
do if (eprin > nfac).
compute eprin=nfac.
end if.
  compute nx2=nx(1,sprin:eprin).
  print pcorr2(1,sprin:eprin)/
     title='偏相関係数'/cname=nx2
     /rname=ng/format=f8.3.

compute sprin=eprin+1.
END LOOP IF (nfac <= eprin).

*--------
*重相関.
*--------.
  compute mcorr=evl(l1)/(evl(l1)+1).
  compute mcorr={mcorr,sqrt(mcorr)}.
  print mcorr/title="     相関比   重相関係数"/format=f8.3.
  compute colname={" ","******","第",num2(1,l1),"相関比解","  ******","part 2"}.
  print colname/format=A8.
  compute sk=0.
  compute rowname={num2}.
  loop i=1 to #nvar.
    compute ek=sk+cmaxs(1,i).
    compute sk=sk+1.
    compute coeffic=t(ncat(1,sk:ek)).
    compute m=dcoef(l1,sk:ek).
    compute coeffic={coeffic,t(m)}.
    compute rlabel=num2.
    compute rlabel(1,1)=nx(1,i).
    compute clabel={"頻度","数量"}.
    print coeffic/format=f8.3/rnames=rlabel/cnames=clabel.  
    compute range=rmax(m)-rmin(m).
    compute r={range, pcorr((#nvar+1),i),corr(#nvar+1,i)}.
    print r/title="               範囲    偏相関係数  相関係数"
             /format=f10.3/rlabels="    --->".
    compute sk=ek.
  end loop. 
end loop.
!if (!ascore2 <> 'NO') !then.
compute vlabel={"h2s1","h2s2","h2s3","h2s4","h2s5","h2s6","h2s7","h2s8","h2s9","10","h2s11","h2s12","h2s13","h2s14","h2s15","h2s16"}.
save dscore/outfile="h2_score.sav"/names=vlabel.
!ifend.
end matrix.
!if (!ascore2 <> 'NO') !then .
match files file="h2_score.sav"/file=*/map.
exec.
erase file="h2_score.sav".
exec.
!ifend.
!if (!dummy2 <> 'NO') !then.
match files file="h2_x.sav"/file=*/map.
exec.
erase file='h2_x.sav'.
exec.
!ifend.
restore.
!enddefine.
restore.
*hayasi2.sps end ====================================================.

*実験用プログラム----------------------------------------------------.
TITLE 数量化２類の検討.
DATA LIST /NUM A1 A2 A3 A4 A5 A6 A7 1-16.
VARIABLE LABELS
 A1 '収入'
 A2 '職業'
 A3 '新しもの好き'
 A4 '広告接触'
 A5 'マスコミ接触'
 A6 '嗜好品購入額（１週間）'
 A7 '小遣い（１週間）'
.
VALUE LABELS
 A1 1 '低所得' 2 '普通' 3 '高所得'/
 A2 1 '事務・技術' 2 '管理・自営' 3 '労務' 4 'その他'/
 A3 1 '好き' 2 '嫌い'/
 A4 1 '見る' 2 '見ない'/
 A5 1 '新聞型' 2 'ＴＶ型' 3 'その他'
.
BEGIN DATA
 1 1 1 1 2 1 210
 2 1 3 2 1 2 1 8
 3 2 1 1 1 1 518
 4 3 4 2 2 2 630
 5 1 4 1 1 1 315
 6 2 3 2 2 2 212
 7 2 1 1 1 1 622
 8 1 3 1 2 2 312
 9 3 2 2 2 3 525
10 1 1 2 1 2 2 8
11 3 1 1 1 31020
12 2 3 1 2 1 410
13 2 1 2 2 2 314
14 1 3 2 2 1 1 6
15 3 2 1 1 3 928
16 2 4 1 1 310 8
17 1 2 2 1 1 310
18 2 1 2 1 3 425
19 3 1 1 2 1 628
20 2 2 1 1 2 628
END DATA.
* (1)基準変数のカテゴリ数が２つの場合----------------------------.
hayasi2 ivar=A1 a3 a2 A5/cvar=A4/pscore=yes/missing=NO.
* (2)基準変数のカテゴリ数が３つの場合----------------------------.
hayasi2 ivar=A1 a3 a2 A4/cvar=A5/ascore=yes/dummy=yes /missing=NO.
* (3)判別分析　(2)の結果をふまえての判別分析---------------------.
discriminant gropus=a5(1,3)/var=col1 to col7/statistics=all/plot=all.
*-------------------------------------------------------------.