c  SEQUENTIAL

c   This program computes a variety of basic sequential analysis statistics.
c   The data are assumed to be a series of integer codes with values ranging
c   from "1" to what ever value the user specifies in the "NCODES"
c   computation at the start of the program.


c  When consecutive codes may not repeat, the expected frequencies
c  and adjusted residuals will be slightly inaccurate for this
c  Fortran 77 program.  This is because iterative proportional fitting
c  and corrections to the computations of adjusted residuals have not
c  been implemented in this program.  See the SAS or SPSS
c  versions of the program for the correct coefficients.


c Enter the appropriate values in the PARAMETER statement for the following:
c NCODES  = number of possible code values
c NEVENTS = number of events, i.e., total number of events in the sequence
c LAG = the lag number for the analyses
c ADJACENT = Enter "1" if adjacent events can be coded the same,
c            enter "0" if not
c TAILED  = Enter "1" for one-tailed tests; enter "2" for two-tailed tests
c PERMTEST = Enter "1" if you want permutation tests of significance,
c            Enter "0" if you don't
c NPERMS  = number of desired permutations per block
c NBLOCKS = number of desired blocks of permutations
c CONFID  = Enter "95" for 95% confidence intervals;
c           "99" for 99% confidence intervals


       parameter ( ncodes   = 6
&,                 nevents  = 122
&,                 lag      = 1
&,                 adjacent = 1
&,                 tailed   = 1
&,                 permtest = 1
&,                 nperms   = 1000
&,                 nblocks  = 10
&,                 confid   = 95      )


      parameter (ncells=ncodes*ncodes)
      integer data (nevents), freqs (ncodes,ncodes), limit,
& results (nperms,ncells), d, kay, signs (ncodes,ncodes),df,
& datap(nevents), datap2(nevents+2), freqsp (ncodes,ncodes),
& obs2 (ncells), signs2 (ncells), coltots (ncodes), perm, block

      real counter, sum, rowtots (ncodes), ntrans,
& pcols (ncodes), p (ncodes,ncodes), prows (ncodes),
& et (ncodes,ncodes), expfreq (ncodes,ncodes),
& zgal (ncodes,ncodes), zadjres (ncodes,ncodes),
& yulesq (ncodes,ncodes), var (ncodes,ncodes),
& min (ncodes,ncodes), kappa (ncodes,ncodes),
& zkappa (ncodes,ncodes), n, nr (ncodes), pnr (ncodes),
& lrx2t, sigs (nblocks,ncells), means (ncells),
& meansigs (ncodes,ncodes), conhi (ncells), conlo (ncells),
& confidhi (ncodes,ncodes), confidlo (ncodes,ncodes), cssq,
& semeans (ncells), nb, obs222 (ncodes,ncodes), obs22 (ncells)


c The program conducts the analyses on a vector
c (a one-dimensional array) of codes.  The name of the
c vector/array must be "data".  The sequences of codes
c can be entered using the DATA or READ commands.


c Entering data using the DATA statement

      data data/3,5,3,4,4,6,3,4,4,1,6,3,6,6,1,6,2,4,3,4,3,4,
& 2,6,6,3,6,3,6,3,4,3,5,3,4,2,5,3,4,2,6,3,4,2,5,3,6,2,6,2,6,
& 3,3,6,3,5,3,6,6,3,3,6,2,6,2,6,3,3,6,3,6,3,3,6,6,2,6,2,6,2,
& 6,3,3,5,3,6,2,6,6,3,3,5,3,5,3,3,5,3,3,6,6,2,6,2,6,2,6,6,2,
& 3,6,2,6,2,6,1,6,2,6,2,6,1 /


c Use these commands to read data from an external file:
c Enter the name of the input data file here, change "unit #" if necessary
c      open (unit=1, file='data', status='old')
c      do 1, i = 1,nevents
c      read (1, 1001) data (i)
c1001  format (i10)
c1     continue

c Enter the name of the output data file here, change "unit #" if necessary
       open (unit=2, file='z',  status='new')


c transitional frequency matrix
       do 2, c = 1,nevents
       if (c + lag .le. nevents) then
       freqs(data(c),data(c+lag)) = freqs(data(c),data(c+lag))+1
       endif
2      continue


c Initializing matrices
       do 3, i = 1,ncodes
       do 4, j = 1,ncodes
       rowtots(i)   =  0 
       coltots(i)   =  0
       prows(i)     =  0
       pcols(i)     =  0
       pnr(i)       =  0
       p(i,j)       = -9999
       et(i,j)      = -9999
       expfreq(i,j) = -9999
       zgal(i,j)    = -9999
       zadjres(i,j) = -9999
       yulesq(i,j)  = -9999
       var(i,j)     = -9999
       min(i,j)     = -9999
       kappa(i,j)   = -9999
       zkappa(i,j)  = -9999
       signs(i,j)   =  0
       meansigs(i,j) = 1
       confidhi(i,j) = 1
       confidlo(i,j) = 1
4      continue
3      continue

       do 5, i = 1,nblocks
       do 6, j = 1,ncells
       sigs(i,j) = 1
6      continue
5      continue

      ntrans = 0
      do 7, i = 1,ncodes
      do 8, j = 1,ncodes
      rowtots(i) = rowtots(i) + freqs(i,j)
      nr(i)      = nr(i)      + freqs(i,j)
      coltots(j) = coltots(j) + freqs(i,j)
      ntrans = ntrans + freqs(i,j)
8     continue
7     continue

      n = ntrans + 1
      nr(data(n)) = nr(data(n)) + 1
      do 9, i = 1,ncodes
      prows(i) = rowtots(i) / ntrans
      pcols(i) = coltots(i) / ntrans
      pnr(i)  = nr(i) / n
9     continue

      lrx2t = 0

      do 10, i = 1,ncodes
      do 11, j = 1,ncodes

      if (adjacent .eq. 0 .and. (ntrans - rowtots(i)) .gt. 0 )   then
      pcols(j) = coltots(j) / (ntrans - rowtots(i) )
      endif
      if (adjacent .eq. 0 .and. (n - nr(i)) .gt. 0)   then
      et(i,j)  = (nr(i) * nr(j)) / (n - nr(i)) 
      endif
      if (adjacent .eq. 0 .and. (ntrans - rowtots(j)) .gt. 0 )   then
      expfreq(i,j)=(rowtots(i)*coltots(j))/(ntrans-rowtots(j))
      endif

      if (adjacent .eq. 1)   then
      et(i,j)  = (nr(i) * nr(j)) / n
      endif
      if (adjacent .eq. 1)   then
      expfreq(i,j) = ( rowtots(i) * coltots(j)) / ntrans 
      endif

c transitional probabilities
      if (rowtots(i) .gt. 0)   then
      p(i,j) = freqs(i,j) / rowtots(i)
      endif

c tablewise LRX2
      if ( freqs(i,j) .gt. 0 .and. expfreq(i,j) .gt. 0) then
      lrx2t = lrx2t + 2 * (freqs(i,j) * log( freqs(i,j)/expfreq(i,j)))
      endif

c Gottman-Allison-Liker z values
      if ( (et(i,j)*(1-pnr(j))*(1-pnr(i))) .gt. 0)   then
      zgal(i,j)=(freqs(i,j)-et(i,j))/sqrt(et(i,j)*(1-pnr(j))*(1-pnr(i)))
      endif

c adjusted residuals (z values) 
      if ( (expfreq(i,j)*(1-pcols(j))*(1-prows(i))) .gt. 0)   then
      zadjres(i,j)=(freqs(i,j)-expfreq(i,j)) /
&     sqrt( expfreq(i,j)*(1-pcols(j)) * (1-prows(i)) )
      endif

c Yule's Q
      a = freqs(i,j)
      b = rowtots(i) - freqs(i,j)
      c = coltots(j) - freqs(i,j)
      d = ntrans - rowtots(i) - coltots(j) + freqs(i,j)
      if ( (a*d + b*c) .gt. 0 )   then
      yulesq(i,j) = ( a*d - b*c ) / (a*d + b*c)
      endif

c kappas & their z values, see Wampold, 1989, & Wampold & Margolin, 1982
      var(i,j)=(nr(i)*nr(j)*(n-nr(j))*(n-nr(i)))/(n**2 *(n-1))
      if (var(i,j) .gt. 0)   then
      zkappa(i,j)  = ( freqs(i,j) - et(i,j) ) / sqrt(var(i,j))
      if ( nr(i) .le. nr(j) ) then
      min(i,j) = nr(i)
      else
      min(i,j) = nr(j)
      endif
      kappa(i,j) = ( freqs(i,j) - et(i,j) ) / ( min(i,j) - et(i,j) )
      if ( kappa(i,j) .lt. 0 ) then
      kappa(i,j) = ( freqs(i,j) - et(i,j) ) /  et(i,j) 
      endif
      endif

c signs
      if ( freqs(i,j) .gt. expfreq(i,j) ) then
      signs(i,j) =  1
      else if ( freqs(i,j) .lt. expfreq(i,j) ) then
      signs(i,j) = -1
      endif

11    continue
10    continue




c Permutation tests of significance
      if (permtest .eq. 1)  then

c reshaping into vectors
      do 100, i = 1,ncodes
      do 101, j = 1,ncodes
      obs2(ncodes*(i-1)+j) = freqs(i,j)
      signs2(ncodes*(i-1)+j) = signs(i,j)
      obs222(i,j) = expfreq(i,j) - (freqs(i,j) - expfreq(i,j))
      obs22(ncodes*(i-1)+j) = obs222(i,j)
101   continue
100   continue


      do 102, block = 1,nblocks

      print *, 'Now computing for block #', block

      do 103, perm  = 1,nperms

c permuting the sequence of codes, alogrithm from Castellan 1992

c when adjacent codes may be the same
      if (adjacent .eq. 1) then
      do 104, i = 1,nevents
      datap(i) = data(i)
104   continue
      do 105, i = 1,(nevents-1)
      kay = int( (nevents - i + 1) * rng(1) + 1 )  + i - 1
      d = datap(i)
      datap(i) = datap(kay)
      datap(kay) = d
105   continue
      endif

c when adjacent codes may NOT be the same
      if (adjacent .eq. 0) then
      do 106, i = 1,nevents+2
      if (i .eq. 1 .or. i .eq. (nevents+2)) then
      datap2(i) = 0
      else
      datap2(i) = data(i-1)
      endif
106   continue
      do 107,  i = 2,(nevents+2-2)
      limit = 10000
      do 108, j = 1,limit
      kay = int(((nevents+2-1) - i + 1) * rng(1) + 1 )  + i - 1
      if ( (datap2(i-1) .ne. datap2(kay)) .and. (datap2(i+1) .ne.
&     datap2(kay)) .and. (datap2(kay-1) .ne. datap2(i)) .and.
&     (datap2(kay+1) .ne. datap2(i)))   go to 109
108   continue
109   d = datap2(i)
      datap2(i) = datap2(kay)
      datap2(kay) = d
107   continue
      do 110, t = 1,nevents
      datap(t) = datap2(t+1)
110   continue
      end if

c transitional frequency matrix for the permuted data
      do 111, i = 1,ncodes
      do 112, j = 1,ncodes
      freqsp (i,j) = 0
112   continue
111   continue

      do 113, c = 1,nevents
      if (c + lag .le. nevents) then
      freqsp(datap(c),datap(c+lag)) = freqsp(datap(c),datap(c+lag))+1
      endif
113   continue


c reshaping new trans freq matrix into a vector & storing in 'results'
      do 114, i = 1,ncodes
      do 115, j = 1,ncodes
      results(perm,((ncodes*(i-1))+j)) = freqsp(i,j)
115   continue
114   continue

103   continue


c sig levels for the current block of permutations

c one-tailed
      if (tailed .eq. 1) then
      do 116, j = 1,ncells
      counter = 0
      do 117, i = 1,nperms
      if      (results(i,j) .ge. obs2(j) .and. signs2(j) .gt. 0) then
      counter = counter + 1
      else if (results(i,j) .le. obs2(j) .and. signs2(j) .lt. 0) then
      counter = counter + 1
      endif
117   continue
      if (signs2(j) .ne. 0) then
      sigs(block,j) =  counter /  nperms
      endif
116   continue
      endif

c two-tailed
      if (tailed .eq. 2) then
      do 118, j = 1,ncells
      counter = 0
      do 119, i = 1,nperms
      if      ( signs2(j) .gt. 0 .and. ((results(i,j) .ge. 
&     obs2 (j)) .or. (results(i,j) .le. obs22(j)))) then
      counter = counter + 1
      else if ( signs2(j) .lt. 0 .and. ((results(i,j) .le.
&     obs2 (j)) .or. (results(i,j) .ge. obs22(j)))) then
      counter = counter + 1
      endif
119   continue
      if (signs2(j) .ne. 0) then
      sigs(block,j) =  counter /  nperms
      endif
118   continue
      endif

102   continue


c  mean significance levels and confidence intervals
      if       (confid .eq. 95 .and. tailed .eq. 1)  then
      z = 1.645
      else if  (confid .eq. 95 .and. tailed .eq. 2)  then
      z = 1.96
      else if  (confid .eq. 99 .and. tailed .eq. 1)  then
      z = 2.326
      else if  (confid .eq. 99 .and. tailed .eq. 2)  then
      z = 2.576
      endif

c mean sig levels
      do 120, j = 1,ncells
      sum = 0
      do 121, i = 1,nblocks
      sum = sum + sigs(i,j)
121   continue
      means(j) = sum / nblocks
120   continue

c standard deviations
      if (nblocks .gt. 1) then
      do 122, j = 1,ncells
      cssq = 0
      do 123, i = 1,nblocks
      cssq = cssq + ((sigs(i,j) - means(j))**2)
123   continue
      nb = nblocks
      semeans(j) = sqrt( (cssq / (nb-1)) )  / sqrt(nb)
122   continue

c confidence intervals
      do 124, j = 1,ncells
      conhi(j) = means(j) + z * semeans(j)
      conlo(j) = means(j) - z * semeans(j)
124   continue
      endif

c reshaping the vectors of means & confid ints into 2-d matrices
      do 125, i = 1,ncodes
      do 126, j = 1,ncodes
      meansigs(i,j) = means(ncodes*(i-1)+j)
      if (nblocks .gt. 1) then
      confidhi(i,j) = conhi(ncodes*(i-1)+j)
      confidlo(i,j) = conlo(ncodes*(i-1)+j)
      endif
126   continue
125   continue

      endif



      write (2,6001)
6001  format (/'Output from Program SEQUENTIAL'/)

      write (2,6002) (ntrans)
6002  format (/'Number of Transitions = ' F10.0)

      write (2,6003) (rowtots(i),i=1,ncodes)
6003  format (/'Row Totals'/6F8.0)

      write (2,6004) (coltots(i),i=1,ncodes)
6004  format (/'Collumn Totals'/6i7)

      if (adjacent .eq. 1) then
      df = (ncodes - 1)**2
      else
      df = (ncodes - 1)**2 - ncodes
      end if
      write (2,6005) lrx2t, df
6005  format (/'Tablewise Likelihood Ratio (Chi-Square)       df'
&     /8x,F10.3,23x,i7)

      write (2,6006)
6006  format (/'"Exp." = Expected Bidirectional Frequencies')

      write (2,6007)
6007  format (/'"tprob." = Transitional Probabilities')

      write (2,6008)
6008  format (/'"zgal" = Gottman-Allison-Liker z values')

      write (2,6009)
6009  format
& (//'Given Target  Freq   Exp.   tprob.   zgal  adjres  YulesQ   Kappa  
& zkappa'/)

      do 132, i = 1,ncodes
      do 133, j = 1,ncodes
      write (2,6010) i,j,freqs(i,j),expfreq(i,j),p(i,j),zgal(i,j),
&                    zadjres(i,j),yulesq(i,j),kappa(i,j),zkappa(i,j)
6010  format(2x,i1,5x,i1,2x,i6,2x,F6.3,2x,F6.3,2x,F6.3,2x,F6.3,2x,F6.3,
&            2x,F6.3,2x,F6.3)
133   continue
132   continue


      if (permtest .eq. 1) then
      write (2,6011) (tailed)
6011  format (//'Requested tail (1 or 2) for Significance Tests =' F4.0)

      write (2,6012) (nperms)
6012  format (/'Number of Permutations Per Block =' i8)

      write (2,6013) (nblocks)
6013  format (/'Number of Blocks of Permutations =' i6)

      write (2,6014) (confid)
6014  format (/'Percentage for the Confidence Intervals =' F4.0)

      write (2,6016)
6016  format (//'Given Target  Freq   MeanSig  ConfidHi ConfidLo'/)

      do 134, i = 1,ncodes
      do 135, j = 1,ncodes
      write (2,6017) i,j,freqs(i,j),meansigs(i,j),confidhi(i,j),
&                    confidlo(i,j)
6017  format(2x,i1,5x,i1,2x,i6,4x,F7.5,2x,F7.5,2x,F7.5)
135   continue
134   continue
      endif

      stop
      end

c Random number generator
      function rng (j)
      integer a, seed
      data a, seed, m / 16807, 1953125, 2147483647 /
      seed = mod(seed*a,m)
      rng = abs(real(seed) / real(m))
      return
      end