* PARADOM - 2.

* This program tests for parallel dominance or "asymmetry in
  predictability," which is the difference in predictability
  between i to j, and j to i (e.g., whether B's behavior is more
  predictable from A's previous behavior than vice versa), as
  described by Wampold (1984, 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 signs   = make(ncodes,ncodes,0).
compute case    = make(ncodes,ncodes,-9999).
compute pd      = make(ncodes,ncodes,-9999).
compute et      = make(ncodes,ncodes,-9999).
compute ett     = make(ncodes,ncodes,-9999).
compute vart    = make(ncodes,ncodes,-9999).
compute min     = make(ncodes,ncodes,-9999).
compute kappa   = make(ncodes,ncodes,-9999).
compute zkappa  = make(ncodes,ncodes,-9999).
compute zeqk    = make(ncodes,ncodes,-9999).
compute obs     = make(ncodes,ncodes,0).
compute pzkappa = make(ncodes,ncodes,1).
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.
compute prow = nr / (msum(nr)).

loop #i = 1 to ncodes.
loop #j = 1 to ncodes.
do if (nr(#i,1) > 0 and nr(#j,1) > 0).
compute pd(#i,#j)=( (freqs(#i,#j) / nr(#i,1)) - prow(#j,1) )/prow(#j,1).
do if (nr(#i,1) > 0 and nr(#j,1) > 0).
compute et(#i,#j)  = (nr(#i,1) * nr(#j,1)) / n.
end if.
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.
end if.
end loop.
end loop.

loop #i = 1 to ncodes.
loop #j = 1 to ncodes.
do if (nr(#i,1) > 0 and nr(#j,1) > 0).

* kappas.
do if ( freqs(#i,#j) = et(#i,#j) ).
compute kappa(#i,#j) = 0.
compute case(#i,#j) = 0.
end if.
* Wampold's 1st case.
do if ( freqs(#i,#j) > et(#i,#j) and freqs(#j,#i) ge et(#j,#i) ).
compute kappa(#i,#j)=(freqs(#i,#j)-freqs(#j,#i)) /
                     (min(#i,#j)-(nr(#i,1)*nr(#j,1)/n)).
compute case(#i,#j) = 1.
end if.
* Wampold's 2nd case.
do if ( freqs(#i,#j) < et(#i,#j) and freqs(#j,#i) le et(#j,#i) ).
compute kappa(#i,#j)=(freqs(#i,#j)-freqs(#j,#i)) /
                     (-1*(nr(#i,1)*nr(#j,1)/n)).
compute case(#i,#j) = 2.
end if.
* Wampold's 3rd case.
do if ( freqs(#i,#j) > et(#i,#j) and freqs(#j,#i) le et(#j,#i) ).
compute kappa(#i,#j)=((freqs(#i,#j)+freqs(#j,#i))-2*(nr(#i,1)*nr(#j,1)/n))/
                     (min(#i,#j) - (nr(#i,1)*nr(#j,1)/n) ).
do if ( kappa(#i,#j) < 0 ).
compute kappa(#i,#j)=((freqs(#i,#j)+freqs(#j,#i))-2*(nr(#i,1)*nr(#j,1)/n))/
                     (nr(#i,1)*nr(#j,1)/n ).
end if.
compute case(#i,#j) = 3.
end if.
* Wampold's 4rth case.
do if ( freqs(#i,#j) < et(#i,#j) and freqs(#j,#i) ge et(#j,#i) ).
compute kappa(#i,#j)=((freqs(#i,#j)+freqs(#j,#i))-2*(nr(#i,1)*nr(#j,1)/n))/
                     (-1* (nr(#i,1)*nr(#j,1)/n) ).
do if (kappa(#i,#j) < 0).
compute kappa(#i,#j)=((freqs(#i,#j)+freqs(#j,#i))-2*(nr(#i,1)*nr(#j,1)/n))/
                     ((nr(#i,1)*nr(#j,1)/n) - min(#i,#j)).
end if.
compute case(#i,#j) = 4.
end if.

* observed frequency, expected frequency, variance, & z.
* same direction.
do if ((pd(#i,#j)>=0 and pd(#j,#i)>=0) or (pd(#i,#j)<=0 and pd(#j,#i)<=0)).
compute obs(#i,#j) = freqs(#i,#j) - freqs(#j,#i) .
compute ett(#i,#j) = 0 .
compute vart(#i,#j)=(2*nr(#i,1)*nr(#j,1)*(n-nr(#i,1)-nr(#j,1)+1))/(n*(n-1)).
compute zkappa(#i,#j)  = obs(#i,#j) / sqrt(vart(#i,#j)).
compute zeqk(#i,#j) = zkappa(#i,#j).
* different directions.
else if ((pd(#i,#j)<=0 and pd(#j,#i)>=0) or 
         (pd(#i,#j)>=0 and pd(#j,#i)<=0)).
compute obs(#i,#j) = freqs(#i,#j) + freqs(#j,#i).
compute ett(#i,#j) = 2*nr(#i,1)*nr(#j,1)/n .
compute vart(#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)=( obs(#i,#j) - ett(#i,#j) ) / sqrt(vart(#i,#j)).
do if (((pd(#i,#j) = 0 or pd(#j,#i) = 0)) and zeqk(#i,#j) < zkappa(#i,#j)).
compute zkappa(#i,#j) = zeqk(#i,#j).
end if.
end if.
compute pzkappa(#i,#j) = (1 - cdfnorm(abs(zkappa(#i,#j)))) * tailed.
end if.

* signs.
do if   ( kappa(#i,#j) > 0 and case(#i,#j) > 0 ).
compute signs(#i,#j) =  1.
else if ( kappa(#i,#j) < 0 and case(#i,#j) > 0 ).
compute signs(#i,#j) = -1.
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 et  /format "f5.4"
  /title="Expected Values/Frequencies"/cnames=b/rnames=b.
print obs /format "f5.4" 
 /title="Observed Parallel Dominance Frequencies"/cnames=b/rnames=b.
print ett /format "f5.4"
 /title="Expected Parallel Dominance Frequencies"/cnames=b/rnames=b.
print case / 
 /title="Sequential Dominance 'Case' Types (Wampold, 1989)"
 /cnames=b/rnames=b.
print /title="Case 1: i increases j, and j increases i".
print /title="Case 2: i decreases j, and j decreases i".
print /title="Case 3: i increases j, and j decreases i".
print /title="Case 4: i decreases j, and j increases i".
print kappa /format "f5.4"
 /title="Parallel Dominance Kappas"/cnames=b/rnames=b.
print { zkappa &* signs } / format "f6.3"
 /title="z values for the Parallel Dominance 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 Parallel Dominance 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).
compute case2   = reshape(case,1,nrow(freqs)*ncol(freqs)).

print signs.


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.

* parallel dominance frequency matrix for permuted data.
compute obsp = make(ncodes,ncodes,0).
loop #i = 1 to ncodes.
loop #j = 1 to ncodes.
do   if (case(#i,#j) = 1 or case(#i,#j) = 2).
compute obsp(#i,#j) = freqsp(#i,#j) - freqsp(#j,#i) .
else if (case(#i,#j) = 3 or case(#i,#j) = 4).
compute obsp(#i,#j) = freqsp(#i,#j) + freqsp(#j,#i) .
end if.
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   ( case2(1,#j) = 1 or case2(1,#j) = 3 ).
do if   ( signs2(1,#j) > 0 and results(#i,#j) >= obs2(1,#j) ).
compute counter = counter + 1.
else if ( signs2(1,#j) < 0 and results(#i,#j) <= obs2(1,#j) ).
compute counter = counter + 1.
end if.
end if.
do if   ( case2(1,#j) = 2 or case2(1,#j) = 4 ).
do if   ( signs2(1,#j) > 0 and results(#i,#j) <= obs2(1,#j) ).
compute counter = counter + 1.
else if ( signs2(1,#j) < 0 and results(#i,#j) >= obs2(1,#j) ).
compute counter = counter + 1.
end if.
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   ( case2(1,#j) = 1 or case2(1,#j) = 3 ).
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 if.
do if   ( case2(1,#j) = 2 or case2(1,#j) = 4 ).
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 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 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.