c TWOCELLS c This program simultaneously tests the unidirectional dependence c of i to j, and the unidirectional dependence of k to L, an c additive pattern described by Wampold and Margolin (1982) c and Wampold (1989, 1992). 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 Enter the appropriate values in the PARAMETER statement for the following: c Enter the appropriate code values for i, j, k, & L for your desired test; 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 ( i = 1 &, j = 6 &, k = 5 &, L = 2 &, 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), & results (nperms), d, kay, sign, limit, & datap(nevents), datap2(nevents+2), freqsp (ncodes,ncodes), & obsp, coltots (ncodes), perm, block real counter, sum, rowtots (ncodes), ntrans, & var, kappa, zkappa, n, nr (ncodes), & sigs (nblocks), means, obs2, obs22, & confidhi, confidlo, cssq, semeans, nb, ett 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, ii = 1,ncodes rowtots(ii) = 0 nr(ii) = 0 coltots(ii) = 0 3 continue ntrans = 0 do 7, ii = 1,ncodes do 8, jj = 1,ncodes rowtots(ii) = rowtots(ii) + freqs(ii,jj) nr(ii) = nr(ii) + freqs(ii,jj) coltots(jj) = coltots(jj) + freqs(ii,jj) ntrans = ntrans + freqs(ii,jj) 8 continue 7 continue n = ntrans + 1 nr(data(n)) = nr(data(n)) + 1 ett = (nr(i)*nr(j) + nr(k)*nr(L) ) /n var=(nr(i)*nr(j)*(n-nr(i))*(n-nr(j))+ & nr(k)*nr(L)*(n-nr(k))*(n-nr(L))+ & 2*nr(i)*nr(j)*nr(k)*nr(L)) / (n**2 *(n-1)) zkappa = ( (freqs(i,j)+freqs(k,L)) - ett ) / sqrt(var) c pzkappa = (1 - cdfnorm(abs(zkappa))) * tailed if (nr(i) .le. nr(j)) then minij = nr(i) else minij = nr(j) endif if (nr(k) .le. nr(L)) then minkL = nr(k) else minkL = nr(L) endif kappa=((freqs(i,j)+freqs(k,L))-(nr(i)*nr(j)+nr(k)*nr(L))/n) & /(minij+minkL-((nr(i)*nr(j)+nr(k)*nr(L))/n)) if ( kappa .lt. 0) then kappa=((freqs(i,j)+freqs(k,L))-(nr(i)*nr(j)+nr(k)*nr(L))/n)/ & (((nr(i)*nr(j)+nr(k)*nr(L))/n)) endif c Permutation tests of significance if (permtest .eq. 1) then obs2 = freqs(i,j)+freqs(k,L) obs22 = ett-(obs2-ett) if (kappa .gt. 0) then sign = 1 else if (kappa .lt. 0) then sign = -1 else sign = 0 endif do 102, block = 1,nblocks print *, 'Now computing for block #', block do 103, perm = 1,nperms c permuting the sequences; alogrithm from Castellan 1992 c when adjacent codes may be the same if (adjacent .eq. 1) then do 104, ii = 1,nevents datap(ii) = data(ii) 104 continue do 105, ii = 1,(nevents-1) kay = int( (nevents - ii + 1) * rng(1) + 1 ) + ii - 1 d = datap(ii) datap(ii) = datap(kay) datap(kay) = d 105 continue endif c when adjacent codes may NOT be the same if (adjacent .eq. 0) then do 106, ii = 1,nevents+2 if (ii .eq. 1 .or. ii .eq. (nevents+2)) then datap2(ii) = 0 else datap2(ii) = data(ii-1) endif 106 continue do 107, ii = 2,(nevents+2-2) limit = 10000 do 108, jj = 1,limit kay = int(((nevents+2-1) - ii + 1) * rng(1) + 1 ) + ii - 1 if ( (datap2(ii-1) .ne. datap2(kay)) .and. (datap2(ii+1) .ne. & datap2(kay)) .and. (datap2(kay-1) .ne. datap2(ii)) .and. & (datap2(kay+1) .ne. datap2(ii))) go to 109 108 continue 109 d = datap2(ii) datap2(ii) = 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, ii = 1,ncodes do 112, jj = 1,ncodes freqsp (ii,jj) = 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 two-cell frequency for permuted data obsp = freqsp(i,j)+freqsp(k,L) results(perm) = obsp 103 continue c sig levels for the current block of permutations c one-tailed if (tailed .eq. 1) then counter = 0 do 121, ii = 1,nperms if (results(ii) .ge. obs2 .and. sign .gt. 0) then counter = counter + 1 else if (results(ii) .le. obs2 .and. sign .lt. 0) then counter = counter + 1 endif 121 continue if (sign .ne. 0) then sigs(block) = counter / nperms endif endif c two-tailed if (tailed .eq. 2) then counter = 0 do 124, ii = 1,nperms if ( sign .gt. 0 .and. ((results(ii) .ge. & obs2) .or. (results(ii) .le. obs22))) then counter = counter + 1 else if ( sign .lt. 0 .and. ((results(ii) .le. & obs2) .or. (results(ii) .ge. obs22))) then counter = counter + 1 endif 124 continue if (sign .ne. 0) then sigs(block) = counter / nperms endif 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 sum = 0 do 126, ii = 1,nblocks sum = sum + sigs(ii) 126 continue means = sum / nblocks c standard deviations if (nblocks .gt. 1) then cssq = 0 do 128, ii = 1,nblocks cssq = cssq + ((sigs(ii) - means)**2) 128 continue nb = nblocks semeans = sqrt( (cssq / (nb-1)) ) / sqrt(nb) c confidence intervals confidhi = means + z * semeans confidlo = means - z * semeans endif endif write (2,6001) 6001 format (/'Output from Program TWOCELLS'/) write (2,6002) 6002 format (//'Transitional Frequencies:') write (2,6003) 6003 format (/'Given Target Freq '/) do 132, ii = 1,ncodes do 133, jj = 1,ncodes write (2,6004) ii,jj,freqs(ii,jj) 6004 format(2x,i1,5x,i1,2x,i6) 133 continue 132 continue write (2,6005) (ntrans) 6005 format (/'Number of Transitions =' F10.0) write (2,6006) (rowtots(ii),ii=1,ncodes) 6006 format (/'Row Totals'/6F8.0) write (2,6007) (coltots(ii),ii=1,ncodes) 6007 format (/'Collumn Totals'/6i7) write (2,6008) 6008 format (/ & 'Simultaneous Two-Cell Test for the following values of i, j, k, & L') write (2,6009) i, j, k, L 6009 format (/' Code i Code j Code k Code L'/(4i9)) write (2,6010) ett, obs2, kappa, zkappa 6010 format (/' Expected Observed kappa z'/(4F10.3)) 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) means 6014 format (/'Mean Significance Level =' F10.6) if (nblocks .gt. 1) then write (2,6015) (confid) 6015 format (/'Percentage for the Confidence Interval =' f4.0) write (2,6016) confidhi 6016 format (/'High End of the Confidence Interval =' F10.6) write (2,6017) confidlo 6017 format (/'Low End of the Confidence Interval =' F10.6) endif endif stop end c Random number generator function rng (j) integer a, seed data a, seed, m / 16807, 1953126, 2147483647 / seed = mod(seed*a,m) rng = abs(real(seed) / real(m)) return end