*  BIDIRECTIONAL - 2.

* 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.

* Would you like to save any matrix output data for use in 
  SPSS or elsewhere? e.g., Would you like to save the transitional
  frequency matrix? the expected frequencies? kappas? z-values?
  If "yes," see the commands at the end of this file.

* Enter/read the data into a matrix with the name "alldata".  The
  program provides results for each session/dyad/group in your data 
  set.  The very first code for each session/dyad/group must
  be > 999.  The program interprets codes > 999 as session/dyad/group
  numbers, and breaks the codes down into dyads/groups on this 
  basis.  The very last code in that you enter, at the very end of 
  the data codes, must be a number > 999 (to signal the end of the 
  data file). Use the following example as a model. It involves 
  data for 3 sessions/dyads/groups; the sessions/dyad/group numbers 
  are 1000, 2000, and 3000; and the last code in the data 
  matrix is 9999999.


compute alldata = {
 1000;
 2;5;2;4;2;5;2;4;2;4;2;2;5;2;4;2;5;3;4;3;4;1;4;4;1;
 2;5;2;1;5;4;2;5;1;4;4;2;3;6;6;2;2;6;5;5;5;6;3;6;3;
 3;6;3;2;5;3;5;2;6;2;2;2;5;2;4;4;5;2;5;2;5;5;4;2;5;
 5;5;2;5;2;5;2;5;2;2;5;2;5;2;2;5;5;2;2;2;5;2;5;2;5;
 2;5;2;2;5;2;5;3;5;2;5;5;2;2;5;2;2;5;2;6;2;5;4;2;5;
 2;4;4;2;1;2;
 2000;
 5;2;5;2;5;5;5;3;5;2;5;5;2;5;5;2;5;2;5;
 2;4;2;4;2;5;4;2;2;5;2;5;5;5;2;5;2;5;5;2;5;2;5;5;2;
 1;2;5;1;5;1;5;5;4;2;2;2;3;6;3;6;3;6;3;6;3;6;3;3;6;
 6;6;5;1;2;5;2;5;5;2;4;5;5;2;5;5;2;5;2;5;2;5;1;5;4;
 2;2;5;2;5;2;5;4;1;4;4;2;4;4;2;4;2;3;5;4;2;2;5;2;6;
 1;4;1;4;5;5;5;4;5;2;4;5;
 3000;
 5;2;5;2;2;2;5;2;2;5;2;2;4;
 2;2;2;6;2;4;3;3;3;3;2;3;6;3;5;3;3;5;3;2;2;2;2;6;3;
 6;3;6;2;6;1;6;3;3;3;3;5;3;5;2;5;3;3;6;3;1;4;1;5;1;
 5;6;3;3;6;3;6;3;6;3;6;3;6;3;3;6;3;6;3;2;5;2;5;2;5;
 2;5;2;2;2;3;5;1;6;1;6;6;6;1;6;6;1;6;1;6;2;5;1;6;6;
 6;2;2;6;6;6;2;6;2;6;2;5;5;5;5;5;5;6;
 9999999 }.




compute dyadnum   = alldata(1,1).
compute data = { 0 }.
compute out = make(1,(ncodes*ncodes+1),0).

* setting up the data for each group/dyad.
loop #dydloop = 1 to nrow(alldata).
do if ( #dydloop > 1 and alldata(#dydloop,1) < 1000 ).
compute data = { data ; alldata(#dydloop,1) }.
end if.
do if ( #dydloop > 1 and alldata(#dydloop,1) ge 1000 ).
compute data = data(2:nrow(data), 1).



* start of group/dyad analyses.



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 /space=3.
print tailed /title="Requested 'tail' (1 or 2) for Significance Tests =".
print dyadnum  /title="Session/Group/Segment Number:".
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 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.


* Would you like to save any matrix output data for use in SPSS
  or elsewhere? e.g., Would you like to save the transitional 
  frequency matrix? the expected frequencies? kappas? z-values?
  If "yes," find the name of the matrix that you would like to
  save (e.g., freqs, ett, kappa, zkappa) and insert this name into 
  the following COMPUTE command.  The command is currently set 
  to save the kappa output matrix.  Substitute an alternative
  matrix name of your choice.

compute matrxout = kappa.


compute out = { out; dyadnum, (reshape(matrxout,1, ncodes*ncodes )) }.
compute data    = {0}.
compute dyadnum   = alldata(#dydloop,1).
end if.
end loop.
compute out = out(2:nrow(out),:).


* The following SAVE command will save the above-specified matrix
  output data in an SPSS system file.  You must provide legal 
  file and directory names for your system 
  (e.g., 'C:\output.sav'), and then "activate" the SAVE 
  command by removing the "*" that currently makes the command 
  a mere comment.

* save out / outfile ='Hard Disk:output.sav'.

print out.


end matrix.


* The requested output matrix has been reshaped into one long row of data.

* get file='Hard Disk:output.sav' / rename (col1=dyadnum) / map .

* descriptives var= all.