*  BIDIRECTIONAL.

* This program tests the bidirectional dependence of 
  behaviors i to j, and j to i, an additive sequential pattern
  described by Wampold and Margolin (1982) and 
  Wampold (1989, 1992).  The data are assumed to be a series
  of integer codes with values ranging from "1" to what ever value the
  user specifies in the "ncodes" computation at the start of the program.
  


*  Set 'seed' to a non-negative integer; 1953125 is good.
SET MXLOOPS=3000000 length=none printback=none width=80  seed = 1953125.
matrix.

* Enter the number of possible code values here.
compute ncodes = 6.

* Enter labels for the codes (5 characters maximum), if desired, here.
compute labels={"Code 1","Code 2","Code 3","Code 4","Code 5","Code 6",
 "Code 7","Code 8","Code 9","Code 10","Code 11","Code 12","Code 13",
 "Code 14","Code 15"}.

* Enter the lag number for the analyses here.
compute lag = 1.

* Can adjacent events be coded the same?  Enter "1" if yes,  enter "0" if no.
compute adjacent = 1.

* Enter "1" for one-tailed tests; enter "2" for two-tailed tests.
compute tailed = 1.

* Do you want to run permutation tests of significance?  "1"=yes, "0"=no;
  Warning: these computations are time-consuming.
compute permtest = 0.

* Enter the number of desired permutations per block here.
compute nperms = 10.

* Enter the number of desired blocks of permutations here.
compute nblocks = 3.

* Confidence intervals for the permutation tests: 
  Enter "95" for 95% confidence intervals;
  enter "99" for 99% confidence intervals.
compute confid = 95.

* Enter/read the data into a matrix with the name "data",
  as in the following example.

compute 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  }.


compute freqs = make(ncodes,ncodes,0).
loop #c = 1 to nrow(data).
do if ( #c + lag le nrow(data) ).
compute freqs((data(#c,1)),(data((#c+lag),1))) =
        freqs((data(#c,1)),(data((#c+lag),1)))  + 1.
end if.
end loop.

* initializing.
compute ett     = make(ncodes,ncodes,-9999).
compute var     = make(ncodes,ncodes,-9999).
compute min     = make(ncodes,ncodes,-9999).
compute kappa   = make(ncodes,ncodes,-9999).
compute zkappa  = make(ncodes,ncodes,-9999).
compute pzkappa = make(ncodes,ncodes,1).
compute signs   = make(ncodes,ncodes,0).
compute obs     = make(ncodes,ncodes,0).
compute rowtots = rsum(freqs).
compute coltots = csum(freqs).
compute n = nrow(data).
compute nr  = rowtots.
compute nr(data(n,1)) = nr(data(n,1)) + 1.

loop #i = 1 to ncodes.
loop #j = 1 to ncodes.
do if (nr(#i,1) > 0 and nr(#j,1) > 0).
compute obs(#i,#j) = freqs(#i,#j)+freqs(#j,#i) .
compute ett(#i,#j)  = 2*nr(#i,1)*nr(#j,1)/n .
compute var(#i,#j) = 2*nr(#i,1)*nr(#j,1)*( nr(#i,1)*nr(#j,1) +
                     (n-nr(#i,1))*(n-nr(#j,1))-n )/(n**2 *(n-1)).
compute zkappa(#i,#j) = ((freqs(#i,#j)+freqs(#j,#i))-ett(#i,#j))
                        /sqrt(var(#i,#j)).
compute pzkappa(#i,#j) = (1 - cdfnorm(abs(zkappa(#i,#j)))) * tailed.
do if ( nr(#i,1) le nr(#j,1) ).
compute min(#i,#j) = nr(#i,1).
else.
compute min(#i,#j) = nr(#j,1).
end if.
compute kappa(#i,#j)=((freqs(#i,#j)+freqs(#j,#i))-ett(#i,#j))
                 / (2 * min(#i,#j) - ett(#i,#j)).
do if ( kappa(#i,#j) < 0 ).
compute kappa(#i,#j)=((freqs(#i,#j)+freqs(#j,#i))-ett(#i,#j)) /
                     ( ett(#i,#j) ).
end if.
do if (nr(#i,1) = nr(#j,1)).
compute kappa(#i,#j)=((freqs(#i,#j)+freqs(#j,#i))-ett(#i,#j))
                 / ((2 * nr(#j,1) -1) - ett(#i,#j)).
end if.

* signs.
do if   ( kappa(#i,#j) > 0 ).
compute signs(#i,#j) =  1.
else if ( kappa(#i,#j) < 0 ).
compute signs(#i,#j) = -1.
end if.

end if.
end loop.
end loop.

compute b = labels(1,1:ncodes).
compute bb = { b,"Totals"}.
print {freqs, rowtots; coltots, msum(rowtots) }  
 /title="Cell Frequencies, Row & Collumn Totals, & N"
 /cnames=bb/rnames=bb.
print obs /format "f5.4" 
 /title="Observed Bidirectional Frequencies"/cnames=b/rnames=b.
print ett /format "f5.4"
 /title="Expected Bidirectional Frequencies"/cnames=b/rnames=b.
print kappa   /format "f5.4" 
 /title="Bidirectional Kappas"/cnames=b/rnames=b.
print zkappa  /format "f6.3" 
 /title="z values for the Bidirectional Kappas"/cnames=b/rnames=b.
print tailed /title="Requested 'tail' (1 or 2) for Significance Tests =".
print pzkappa /format "f5.4" /title=
 "Significance Levels for the Bidirectional Kappas"/cnames=b/rnames=b.




* Permutation tests of significance.
do if (permtest = 1).

compute obs2    = reshape(obs,1,nrow(freqs)*ncol(freqs)).
compute obs22   = reshape((ett-(obs-ett)),1,nrow(freqs)*ncol(freqs)).
compute signs2  = reshape(signs,1,nrow(freqs)*ncol(freqs)).
compute sigs    = make(nblocks,nrow(freqs)*ncol(freqs),1).

loop #block = 1 to nblocks.
print #block /title="Currently computing for block #:".
print /space=1.

compute results = make(nperms,nrow(freqs)*ncol(freqs),-9999).

loop #perm = 1 to nperms.

* permuting the sequences; alogrithm from Castellan 1992.

* when adjacent codes may be the same.
compute datap = data.
do if (adjacent = 1).
loop #i = 1 to (nrow(datap) -1).
compute kay=trunc( (nrow(datap) - #i + 1) * uniform(1,1) + 1 ) + #i - 1.
compute d = datap(#i,1).
compute datap(#i,1) = datap(kay,1).
compute datap( kay,1) = d.
end loop.
end if.

* when adjacent codes may NOT be the same.
do if (adjacent = 0).
compute datap = { 0; data; 0 }.
loop #i = 2 to (nrow(datap) - 2).
compute limit = 10000.
loop #j = 1 to limit.
compute kay=trunc(((nrow(datap)-1)-#i+1) * uniform(1,1) + 1 ) + #i - 1.
do if 
    ((datap(#i-1,1)  ne datap(kay,1)) and (datap(#i+1,1)  ne datap(kay,1)) 
 and (datap(kay-1,1) ne datap(#i,1))  and (datap(kay+1,1) ne datap(#i,1))).
break.
end if.
end loop.
compute d = datap(#i,1).
compute datap(#i,1) = datap(kay,1).
compute datap( kay,1) = d.
end loop.
compute datap = datap(2:(nrow(datap)-1),:).
end if.


* transitional frequency matrix for permuted data.
compute freqsp = make(ncodes,ncodes,0).
loop #c = 1 to nrow(datap).
do if ( #c + lag le nrow(datap) ).
compute freqsp((datap(#c,1)),(datap((#c+lag),1))) =
        freqsp((datap(#c,1)),(datap((#c+lag),1)))  + 1.
end if.
end loop.

* bidirectional frequency matrix for permuted data.
compute obsp = make(ncodes,ncodes,0).
loop #i = 1 to ncodes.
loop #j = 1 to ncodes.
compute obsp(#i,#j) = freqsp(#i,#j)+freqsp(#j,#i) .
end loop.
end loop.

compute results(#perm,:) =  reshape(obsp,1,nrow(freqs)*ncol(freqs)).
end loop.


* sig levels for the current block of permutations.

* one-tailed.
do if (tailed = 1).
loop #j = 1 to ncol(results).
compute counter = 0.
loop #i = 1 to nrow(results).
do if   ( results(#i,#j) >= obs2(1,#j) and signs2(1,#j) > 0 ).
compute counter = counter + 1.
else if ( results(#i,#j) <= obs2(1,#j) and signs2(1,#j) < 0 ).
compute counter = counter + 1.
end if.
end loop.
do if (signs2(1,#j) ne 0).
compute sigs(#block,#j) = counter /  nperms.
end if.
end loop.
end if.

* two-tailed.
do if (tailed = 2).
loop #j = 1 to ncol(results).
compute counter = 0.
loop #i = 1 to nrow(results).
do if   ( signs2(1,#j) > 0 and ((results(#i,#j) >= obs2 (1,#j)) or
                                (results(#i,#j) <= obs22(1,#j)))).
compute counter = counter + 1.
else if ( signs2(1,#j) < 0 and ((results(#i,#j) <= obs2 (1,#j)) or
                                (results(#i,#j) >= obs22(1,#j)))).
compute counter = counter + 1.
end if.
end loop.
do if (signs2(1,#j) ne 0).
compute sigs(#block,#j) = counter /  nperms.
end if.
end loop.
end if.

end loop.


* mean significance levels and confidence intervals.
do if    (confid = 95 and tailed = 1).
compute z = 1.645.
else if  (confid = 95 and tailed = 2).
compute z = 1.96.
else if  (confid = 99 and tailed = 1).
compute z = 2.326.
else if  (confid = 99 and tailed = 2).
compute z = 2.576.
end if.

compute meansigs = csum(sigs) / nblocks.
do if (nblocks > 1).
compute semeans  = make(1, ncol(sigs), -9999).
loop #a = 1 to ncol(sigs).
compute semeans(1,#a) = (sqrt(cssq (sigs(:,#a)-meansigs(1,#a)) /
                        (nblocks -1))) / (sqrt(nblocks)).
end loop.
compute confidhi = meansigs + z &* semeans.
compute confidlo = meansigs - z &* semeans.
end if.

print          /title="Permutation Tests of Significance:".
print nperms   /title="Number of permutations per block: ".
print nblocks  /title="Number of blocks of permutations: ".
print {reshape(meansigs,ncodes,ncodes)} /format="f7.5"
 /title="Mean Significance Levels"/cnames=b/rnames=b.
do if (nblocks > 1).
print confid /title="Percentage for the Confidence Intervals:".
print {reshape(confidhi,ncodes,ncodes)} /format="f7.5"
 /title="High Ends of the Confidence Intervals"/cnames=b/rnames=b.
print {reshape(confidlo,ncodes,ncodes)} /format="f7.5"
 /title="Low Ends of the Confidence Intervals"/cnames=b/rnames=b.
end if.

end if.


end matrix.