c NONPARALLEL.DOM c This program tests for nonparallel dominance, a form of asymmetry c in predictability, between i to j, and k to L as described by c Wampold (1984, 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 = 5 &, j = 3 &, k = 6 &, 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, case, 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), prow (ncodes), & sigs (nblocks), means, obs2, np, nrp (ncodes), & confidhi, confidlo, cssq, semeans, nb, ett, & et (ncodes,ncodes), pd (ncodes, ncodes), obs22 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 do 4, jj = 1,ncodes rowtots(ii) = 0 nr(ii) = 0 coltots(ii) = 0 pd(ii,jj) = -9999 et(ii,jj) = -9999 4 continue 3 continue do 5, ii = 1,nblocks sigs(ii) = 1 5 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 do 9, ii = 1,ncodes prow(ii) = nr(ii) / ntrans 9 continue n = ntrans + 1 nr(data(n)) = nr(data(n)) + 1 do 10, ii = 1,ncodes do 11, jj = 1,ncodes if (nr(ii) .gt. 0 .and. nr(jj) .gt. 0 .and. & prow(jj) .gt. 0) then pd(ii,jj) = ( (freqs(ii,jj) / nr(ii)) - prow(jj) ) / prow(jj) if (nr(ii) .gt. 0 .and. nr(jj) .gt. 0) then et(ii,jj) = (nr(ii) * nr(jj)) / n endif endif 11 continue 10 continue 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 c kappa kappa = -9999 case = -9999 if ( freqs(i,j) .eq. et(i,j) ) then kappa = 0 case = 0 endif if ( nr(i) .gt. 0 .and. nr(j) .gt. 0 .and. nr(k) .gt. 0 & .and. nr(L) .gt.0 ) then c Wampold's 1st case if ( freqs(i,j) .gt. et(i,j) .and. freqs(k,L) & .ge. et(k,L) ) then kappa=(nr(k)*nr(L)*freqs(i,j)-nr(i)*nr(j)*freqs(k,L))/ & ( nr(k)*nr(L)*minij - (nr(i)*nr(j)*nr(k)*nr(L)/n)) if ( kappa .lt. 0 ) then kappa=(nr(k)*nr(L)*freqs(i,j)-nr(i)*nr(j)*freqs(k,L))/ & ( nr(i)*nr(j)*minkL - (nr(i)*nr(j)*nr(k)*nr(L)/n)) endif case = 1 endif c Wampold's 2nd case if ( freqs(i,j) .lt. et(i,j) .and. freqs(k,L) & .le. et(k,L) ) then kappa=(nr(k)*nr(L)*freqs(i,j)-nr(i)*nr(j)*freqs(k,L))/ & (-1*(nr(i)*nr(j)*nr(k)*nr(L)/n)) case = 2 endif c Wampold's 3rd case if ( freqs(i,j) .gt. et(i,j) .and. freqs(k,L) & .le. et(k,L) ) then kappa=(nr(k)*nr(L)*freqs(i,j)+nr(i)*nr(j)*freqs(k,L)- & 2*(nr(i)*nr(j)*nr(k)*nr(L)/n) )/ & ( nr(k)*nr(L)*minij - (nr(i)*nr(j)*nr(k)*nr(L)/n)) if ( kappa .lt. 0 ) then kappa=(nr(k)*nr(L)*freqs(i,j)+nr(i)*nr(j)*freqs(k,L)- & 2*(nr(i)*nr(j)*nr(k)*nr(L)/n) )/ (nr(i)*nr(j)*nr(k)*nr(L)/n) endif case = 3 endif c Wampold's 4rth case if ( freqs(i,j) .lt. et(i,j) .and. freqs(k,L) .ge. & et(k,L) ) then kappa=(nr(k)*nr(L)*freqs(i,j)+nr(i)*nr(j)*freqs(k,L)- & 2*(nr(i)*nr(j)*nr(k)*nr(L)/n) )/(-1*(nr(i)*nr(j)*nr(k)*nr(L)/n)) if (kappa .lt. 0) then kappa=(nr(k)*nr(L)*freqs(i,j)+nr(i)*nr(j)*freqs(k,L)- & 2*(nr(i)*nr(j)*nr(k)*nr(L)/n) )/ & ((nr(i)*nr(j)*nr(k)*nr(L)/n) - nr(i)*nr(j)*minkL ) endif case = 4 endif endif c observed frequency, expected frequency, variance, & z zeqk = 9999 obs = -9999 ett = -9999 zkappa = -9999 pzkappa = -9999 if (pd (i,j) .ne. -9999 .and. pd(k,L) .ne. -9999) then c same direction if ( (pd(i,j) .ge. 0 .and. pd(k,L) .ge. 0) .or. & (pd(i,j) .le. 0 .and. pd(k,L) .le. 0) ) then obs=nr(k)*nr(L)*freqs(i,j)-nr(i)*nr(j)*freqs(k,L) ett = 0 var=(nr(i)*nr(j)*nr(k)*nr(L)* & (n*nr(i)*nr(j) + n*nr(k)*nr(L) - & nr(i)*nr(j)*nr(k) - nr(i)*nr(j)*nr(L) - & nr(i)*nr(k)*nr(L) - nr(j)*nr(k)*nr(L) ))/ (n*(n-1)) zkappa = obs / sqrt(var) zeqk = zkappa c different directions else if ( (pd(i,j) .le. 0 .and. pd(k,L) .ge. 0) & .or. (pd(i,j) .ge. 0 .and. pd(k,L) .le. 0) ) then obs=nr(k)*nr(L)*freqs(i,j)+nr(i)*nr(j)*freqs(k,L) ett = 2*nr(i)*nr(j)*nr(k)*nr(L)/n var=(nr(i)*nr(j)*nr(k)*nr(L)* & (4*nr(i)*nr(j)*nr(k)*nr(L) + n**2 *nr(i)*nr(j) + & n**2 *nr(k)*nr(L) - n*nr(i)*nr(j)*nr(k) - n*nr(i)*nr(j)*nr(L)- & n*nr(i)*nr(k)*nr(L) - n*nr(j)*nr(k)*nr(L) )) / (n**2 *(n-1)) zkappa = ( obs - ett ) / sqrt(var) if ( ((pd(i,j) .eq. 0 .or. pd(j,i) .eq. 0)) .and. & zeqk .lt. zkappa ) then zkappa = zeqk endif endif c pzkappa = (1 - cdfnorm(abs(zkappa))) * tailed endif if (kappa .gt. 0) then sign = 1 else if (kappa .lt. 0 .and. kappa .ne. -9999) then sign = -1 else sign = 0 endif c Permutation tests of significance if (permtest .eq. 1 .and. obs .ne. -9999 .and. & ett .ne. -9999) then obs2 = obs obs22 = ett-(obs2-ett) 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 nonparallel dominance frequency for permuted data np = nevents do 1135, ii = 1,ncodes nrp(ii) = 0 1135 continue do 114, ii = 1,ncodes do 115, jj = 1,ncodes nrp(ii) = nrp(ii) + freqsp(ii,jj) 115 continue 114 continue nrp(datap(np)) = nrp(datap(np)) + 1 if (case .eq. 1 .or. case .eq. 2) then obsp=nrp(k)*nrp(L)*freqsp(i,j)-nrp(i)*nrp(j)*freqsp(k,L) else if (case .eq. 3 .or. case .eq. 4) then obsp=nrp(k)*nrp(L)*freqsp(i,j)+nrp(i)*nrp(j)*freqsp(k,L) endif 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 ( case .eq. 1 .or. case .eq. 3 ) then if ( sign .gt. 0 .and. results(ii) .ge. obs2 ) then counter = counter + 1 else if ( sign .lt. 0 .and. results(ii) .le. obs2 ) then counter = counter + 1 endif endif if ( case .eq. 2 .or. case .eq. 4 ) then if ( sign .gt. 0 .and. results(ii) .le. obs2 ) then counter = counter + 1 else if ( sign .lt. 0 .and. results(ii) .ge. obs2 ) then counter = counter + 1 endif 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 ( case .eq. 1 .or. case .eq. 3 ) then 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 endif if ( case .eq. 2 .or. case .eq. 4 ) then if ( sign .gt. 0 .and. ((results(ii) .le. obs2) .or. & (results(ii) .ge. obs22))) then counter = counter + 1 else if ( sign .lt. 0 .and. ((results(ii) .ge. obs2) .or. & (results(ii) .le. obs22))) then counter = counter + 1 endif 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 NONPARADOM'/) write (2,6002) 6002 format (//'Transitional Frequencies & Expected Frequencies') write (2,6003) 6003 format (/'Given Target Frequency Expected'/) do 132, ii = 1,ncodes do 133, jj = 1,ncodes write (2,6004) ii,jj,freqs(ii,jj),et(ii,jj) 6004 format(2x,i1,5x,i1,5x,i6, 3x,F8.3) 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 (// & 'Nonparallel Dominance 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) case 6010 format (/'Sequential Dominance Case Type (Wampold, 1989) = ' i6) write (2,6013) 6013 format( & /'Case 1: i increases j, and j increases i' & /'Case 2: i decreases j, and j decreases i' & /'Case 3: i increases j, and j decreases i' & /'Case 4: i decreases j, and j increases i'/) write (2,6011) ett, obs, kappa, (zkappa*sign) 6011 format (/' Expected Observed kappa z'/(4F10.3)) write (2,6012) 6012 format (/ & 'The above Expected & Observed values are weighted; see Wampold, 1989, p &184'/) if (permtest .eq. 1 .and. obs .ne. -9999 .and. & ett .ne. -9999) then write (2,6014) (tailed) 6014 format (/'Requested tail (1 or 2) for Significance Tests =' F4.0) write (2,6015) (nperms) 6015 format (/'Number of Permutations Per Block =' i8) write (2,6016) (nblocks) 6016 format (/'Number of Blocks of Permutations =' i6) write (2,6017) means 6017 format (/'Mean Significance Level =' F10.6) if (nblocks .gt. 1) then write (2,6018) (confid) 6018 format (/'Percentage for the Confidence Interval =' f4.0) write (2,6019) confidhi 6019 format (/'High End of the Confidence Interval =' F10.6) write (2,6020) confidlo 6020 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, 1953125, 2147483647 / seed = mod(seed*a,m) rng = abs(real(seed) / real(m)) return end