$sub QSOD
  !
  !------------------------------------------------------------------
  !
  !   Author:  M. P. Becker,  Department of Biostatistics
  !                           University of Michigan, U.S.A.
  !   Main macros:
  !        QSOD   Fits multiplicative interaction models to
  !               square contingency tables.
  !        SCALE  Scales scores estimated with QSYM so that they satisfy
  !               the restrictions NU(%1) = %1 and NU(%2) = %2.
  !
  !        For macro QSYM:
  !               Macros required:
  !                 LINP  The LINear Part of the model.
  !
  !               Variables required:
  !                 NU    Initial estimates of the scores
  !                 IND   Indicator variable specifying which scores
  !                       are to be estimated.
  !
  !               Scalar Arguements:
  !                 %R    Dimensions of the table (%r X %r).
  !                 %W    Maximum number of iterations (default 10).
  !
  !               Variables deleted and/or created:
  !                 r_, c_, beta, f_, n2_, n3_, x_, xx_, y_, z_, q_
  !                 DF_, NE_
  !
  !               Scalars used:
  !                 %b, %d, %e, %f, %h, %i, %j, %r, %s,
  !                 %x, %w, %z1, %z2, %z3, %z4
  !-------------------------------------------------------------------
  $macro QSOD! QuasiSymmetry macros for Ordinal Data
  $warn$
  $delete r_ c_ beta f_ n2_ n3_ x_ xx_ y_ z_ q_$
  $calc r_=%gl(%r,%r) : c_=%gl(%r,1) : %h=1 : %d=0$
  $factor r_ %r c_ %r$
  $var %r n2_ x_ f_ n3_$
  $var 1 DF_ NE_ $
  $calc NE_=%cu(IND)                     ! number of scores estimated
  $calc %z1=1 : %z2=%coc : %z3=%if(%gt(%w,0),%w,10)$
  $calc beta = nu(r_) * nu(c_)$
  $print : 'Deviance     df     Iteration' $
  $out $
  $fit #linp+beta$
  $out %z2 $
  $use depr $
  $warn $
  $extract %pe$
  $calc %b=%pe(%pl)$
  $recycle $
  $while %z1 FITQ$
  $cycle $
  $endmac
  !
  $macro FITQ! FIT Quasisymmetric models for ordinal data
  $warn $
  $calc %h=%h+1$
  $calc %j=1 : %s=1$
  $while %s UP$
  $calc n2_ = nu : nu = nu + %eq(ind,1) * x_ : beta = nu(r_) * nu(c_)$
  $out $
  $fit .  $
  $out %z2 $
  $use depr $
  $calc %f = %b $
  $extract %pe $
  $calc %b=%pe(%pl) $
  $calc f_=(n2_-nu)**2 $sort n3_ f_ $calc %d=n3_(%r)$
  $calc %d=%sqrt(%d)$
  $calc %f=%sqrt((%f-%b)**2)$
  $calc %d=%if(%lt(%f,%d),%d,%f)$
  $calc %i=%lt(%d,0.001)$
  $switch %i CONV$
  $calc %e=%ge(%h,%z3)$
  $switch %e MADE$
  $warn $
  $endmac
  !
  $macro UP!   UPdate the scores
  $calc xx_ = (%eq(%j,r_)*nu(c_)*(%yv-%fv))+(%eq(%j,c_)*nu(r_)*(%yv-%fv))$
  $calc xx_ = xx_*%b$
  $calc %x = %cu(xx_)$
  $calc y_ = %ne(r_,c_)*((%eq(%j,r_)*(nu(c_)**2)*%fv)
                        +(%eq(%j,c_)*(nu(r_)**2)*%fv))$
  $calc y_ = y_ * %b**2 $
  $calc z_ = %eq(r_,c_)*(%eq(%j,r_)*(4*(nu(%j)**2)*(%b**2)*%fv))$
  $calc q_ = %eq(r_,c_)*(%eq(%j,r_)*(2*%b*(%yv-%fv)))$
  $calc y_ = y_ + z_ - q_$
  $calc %y = %cu(y_)$
  $calc x_(%j) = %x/%y : %j=%j+1$
  $calc %s=%le(%j,%r)$
  $endmac
  !
  $macro DEPR ! DEviance PRint
  $calc DF_=%df-NE_+(2-%r+NE_)*%lt((%r-NE_),2) ! corrected df $
  $print *r %dv,8,3  *i DF_,7  *i %h,10 $
  $endmac
  !
  $macro MADE!  MAx absolute Diff (change) in Estimated scores
  $print 'ALGORITHM DID NOT CONVERGE IN '*i %z3,4 ' CYCLES OF ITERATIONS'$
  $print 'MAXIMUM ABSOLUTE CHANGE IN ESTIMATED SCORES:'%d$
  $print ' '$
  $calc %z1=0$
  $endmac
  !
  $macro CONV!  check for CONVergence
  $print ; '    SCORES'$
  $look nu$
  $calc %z4=%ne(%df,DF_) $
  $switch %z4 warn $
  $calc beta=nu(r_)*nu(c_) : %z1=0 $fit . $dis e$
  $endmac
  !
  $macro WARN ! WARNing message regarding s.e.'s and DF
  $print ; 'Standard errors and d.f. given below are not valid.' ; $
  $endmac $
  !
  $return !
  !
  $sub scale !
  $macro SCALE ! SCALE scores so that nu(%1)=%1 and nu(%2)=%2
  !
  !      Formal arguements:
  !              %1 & %2   Categories to have scores fixed at their
  !                        number.
  !
  $warn $
  $calc %a=(%2-%1)/(nu(%2)-nu(%1)) : %b=%1 - (%a * nu(%1))$
  $calc nu = %a * nu + %b$
  $use conv $
  $warn $
  $endmac
  !
  $return
  !
  $finish$