'p値を用いたステップワイズ(プロビットモデルの場合)

'その2
'その1として作成した一番大きなp値を持つ変数を見つける
'プログラムに、変数を取り除いたうえで、閾値以下に
'なる変数だけを残すようにWhile loopで拡張しました。


'ここから---------------------------------------
'推定したオブジェクト名を%eqnameとします。
'ここでは例としてプロビットモデルのeq01を利用します。
%eqname="eq01"
'p値の境界値を設定します。0.1とした場合、
'一番悪い(大きな)p値として0.1を設定した事になります。
!cp=0.1

'元の推定式を***orgという名前でコピーしておきます
%newname=%eqname+"org"
copy %eqname %newname
'----------------------------------------------
'このプログラムは説明変数が2個以上で機能します。
'----------------------------------------------
'p値は直接データメンバとして用意されていませんので、
'z値と累積分布関数から計算し、p値をシリーズcheckpに
'コピーします。ただし、定数項のp値は操作の対象外とします。

'ループを開始するため、仮に!maxをゼロとします。
'ループ開始後に!maxには一番大きなp値が入ります。
!max=1

while !max > !cp

!num=0
!noc={%eqname}.@ncoef

'----------------------------------------------
'説明変数の数(定数項を含まない)が2個しかない場合は
'ステップワイズは実行しません(その2で追加)。
if !noc < 3 then 
@uiprompt("Error: Only two independent variables") 
stop
else
endif
'----------------------------------------------
series checkp
checkp=na
for !i=2 to !noc
'z値の符号によって計算式が異なるので、z値は絶対値を
'とります。
!z=@abs({%eqname}.@tstats(!i))
!p=(1-@cnorm(!z))*2
checkp(!i-1)=!p
next
'これで推定式からp値を取り出して、計算用のcheckpという
'シリーズの準備ができました。

'次に一番大きなp値を持つ変数を求めます。

!max=checkp(1)
for !j=1 to !noc-1
if checkp(!j+1) > !max then
	!max=checkp(!j+1)
	!num=!j+2
else
endif
next

'!maxが!cp(pの境界値)よりも小さい場合、以下のループは実行
'しません。
'主にここから下のコードを「その2」で追加しました。
'-------------------------------------------------
if !max < !cp then
else

'変数名(series object)を取得する
%indv=eq01.@varlist
%dep=@wmid(%indv,1,1)

'定数項(2番目)から外す変数の1つ手前までの文字列を取得
%s1=@wmid(%indv,2,!num-1)

'説明変数の個数
!ncoef=eq01.@ncoef
'外す変数の後ろから最後までの文字列を取得
%s2=@wmid(%indv,!num+2,!ncoef-!num)
%ns=%dep+" "+%s1+" "+%s2
%mtd=eq01.@method
'外した後の推定式はeq01を上書きします。
equation eq01.{%mtd} {%ns}
endif

wend
show eq01