$TRAN i w f h o
$c This was prepared for Class 5, Soc. 619, Fall 1993
$C Data from Clogg, Table 5
$C 
$UNITS 20 $DATA FREQ $READ
44 11 38 62
59 41 147 293
23 11 13 27
27 8 16 27
258 57 105 110
$CAL %r=5
$CAL %c=4
$var %r mu $var %c nu $
$calc mu=%GL(%r,1)
$CAL nu=%GL(%c,1)  !initital parameter estimates $

  $macro RC !
  $delete r_ c_ m1_ m2_  f1_ f2_ f3_ n1_ n2_ n3_  df_ beta $
  $warn !
  $calc %z1=%if(%gt(%w,0),%w,10) $
  $calc r_=%gl(%r,%c) : c_=%gl(%c,1) : %h=0 : %d=0  $
  $var 1 df_ $
  $factor r_ %r c_ %c  $
  $calc %z=1 : %z2=%coc  $
  $calc f2_=%yv : n3_=0 $
  $print : 'Deviance     df     Iteration     Change' $
  $while %z frc          $
  $delete df_ f1_ f2_ f3_ m1_ m2_ n1_ n2_ n3_  r_ c_  beta $
  $warn !
  $endmac                !
  !
  $macro frc             !
  $warn                  $
  $out $tran $
  $calc n2_=nu(c_)      $
  $fit #linp+r_.n2_     $
  $var  %pl m1_              $
  $extract %pe $calc m1_=%pe $
  $calc %m=%pl-%r            $
  $calc mu=%gl(%r,1) : mu(mu)=m1_(mu+%m) : m2_=mu(r_) $
  $fit #linp+c_.m2_  $
  $calc %h=%h+1   $
  $var  %pl n1_        $
  $extract %pe $calc n1_=%pe      $
  $calc nu=%gl(%c,1) : %o=%pl-%c :  nu(nu)=n1_(nu+%o) $
  $calc f1_=f2_ : f2_=%fv        $
  $calc f3_=(f2_-f1_)**2 $
  $sort n3_ f3_ $calc %d=n3_(%nu) $
  $calc %d=%sqrt(%d)  $
  $out %z2 $tran i o $
  $use depr $
  $calc %d=%lt(%d,0.001) $
  $switch %d last $
  $calc %e=%ge(%h,%z1)  $
  $switch %e drop   $
  $warn $
  $endmac        !
  !
  $macro depr !
  $calc df_=(%r-2)*(%c-2) $
  $print *r %dv,8,3 *i df_,7 *i %h,10  *r %d,15,3 $
  $endmac
  !
  $macro drop   !
  $print  ;'ALGORITHM DID NOT CONVERGE IN '*i %z1,4 ' CYCLES OF ITERATIONS'  !
  $print ' '$
  $calc %z=0 $
  $endmac      !
  !
  $macro last  !
  $warn $
  $print ' '   !
  $print 'ROW SCORES' !
  $calc %x=%cu(mu)/%r : mu=mu-%x : %x=%cu(mu**2) : mu=mu/%sqrt(%x) !
  $calc %y=%cu(nu)/%c : nu=nu-%y : %y=%cu(nu**2) : nu=nu/%sqrt(%y) !
  $look mu    !
  $print ' '   !
  $print 'COLUMN SCORES' !
  $look nu !
  $calc m2_=mu(r_) : n2_=nu(c_) : %h=0  !
  $print  ;'Standard errors and d.f. given below are not valid. '
  $print 'Correct d.f. = '  *i df_,3$
  $print  ;' ' $
  $calc beta=m2_*n2_ : %z=0 $fit r_+c_+beta $dis e !
  $warn !
  $endmac !
  !
$macro linp r_+c_ $endmacro
$error p
$yvar freq
$c  Independence model
$calc r_=%gl(%r,%c) : c_=%gl(%c,1) 
$factor r_ %r c_ %c $
$fit r_+ c_ 
$disp e$
$use rc $
$stop