!    Author:  M. P. Becker, Department of Biostatistics
  !                           University of Michigan, U.S.A.
  !
  !    RC: Fits multiplicative interaction models to R-by-C
  !        contingency tables.
  !          i.e., log YV(i,j) = LINearPart + BETA * MU(I) * NU(J)
  !
  !    For macro RC:
  !            Macros required:
  !               LINP  The LINear Part of the model.
  !
  !            Variables required:
  !               MU  Initial estimates of the row scores (length %R).
  !               NU  Initial estimates of the column scores (length %C).
  !
  !            Scalar arguments:
  !               %R  Number of rows in the table.
  !               %C  Number of columns in the table.
  !
  !            Variables created and/or deleted:
  !               R_, C_, BETA, F1_, F2_, F3_, N1_, N2_, N3_,
  !               M1_, M2_, DF_
  !
  !            Scalars used:
  !               %c, %d, %e, %h, %o, %r, %w, %x, %y, %z, %z1, %z2
  !
  !       example:
  !         $calc %R = # of rows, %C = # of columns $
  !         $var %R MU $var %C NU $
  !         $calc MU=%gl(%R,1) : NU=%gl(%C,1) ! initial parameter estimates $
  !         $macro LINP R_+C_ $endmac  ! LINearPart of the model
  !         $use RC $
  !
  !
  !--------------------------------------------------------------------
  $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 !
  !
  $return !
  $finish