*  TWOCELLS.

* This program simultaneously tests the unidirectional dependence
  of i to j, and the unidirectional dependence of k to L, an
  additive 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.

* This program simultaneously tests the unidirectional independence
  of i to j, and the unidirectional independence of k to L;
  Enter the code values for your desired test below.
compute i =  1  .
compute j =  6  . 
compute k =  5  .
compute L =  2  .

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

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 ett = (nr(i,1)*nr(j,1) + nr(k,1)*nr(L,1) ) /n .
compute var=(nr(i,1)*nr(j,1)*(n-nr(i,1))*(n-nr(j,1))+
            nr(k,1)*nr(L,1)*(n-nr(k,1))*(n-nr(L,1))+
            2*nr(i,1)*nr(j,1)*nr(k,1)*nr(L,1)) / (n**2 *(n-1)).
compute zkappa  = ( (freqs(i,j)+freqs(k,L)) - ett ) / sqrt(var).
compute pzkappa = (1 - cdfnorm(abs(zkappa))) * tailed.
do if (nr(i,1) le nr(j,1)).
compute minij  = nr(i,1).
else.
compute minij  = nr(j,1).
end if.
do if (nr(k,1) le nr(L,1)).
compute minkL  = nr(k,1).
else.
compute minkL  = nr(L,1).
end if.
compute kappa=((freqs(i,j)+freqs(k,L))-(nr(i,1)*nr(j,1)+nr(k,1)*nr(L,1))/n) 
              /(minij+minkL-((nr(i,1)*nr(j,1)+nr(k,1)*nr(L,1))/n)).
do if ( kappa < 0).
compute kappa=((freqs(i,j)+freqs(k,L))-(nr(i,1)*nr(j,1)+nr(k,1)*nr(L,1))/n)/
              (((nr(i,1)*nr(j,1)+nr(k,1)*nr(L,1))/n)).
end if.

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 { i, j, k, L } /format "f12.4"/title = 
 "Simultaneous Two-Cell Test for the following values of i, j, k, & L"
 /cnames={"Code i =","Code j =", "Code k =", "Code L =" }.
print tailed /title="Requested 'tail' (1 or 2) for Significance Tests =".
print {ett, (freqs(i,j)+freqs(k,L)), kappa, zkappa, pzkappa }
 /format "f12.4"
 /title="Results:"/cnames={"Expected", "Observed", "kappa", "z", "sig."}.




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

compute obs2    = freqs(i,j)+freqs(k,L).
compute obs22   = ett-(obs2-ett).
do if   (kappa > 0).
compute sign =  1.
else if (kappa < 0).
compute sign = -1.
else.
compute sign = 0.
end if.
compute sigs = make(nblocks,1,1).

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

compute results = make(nperms,1,-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.

* two-cell frequency for permuted data.
compute obsp = freqsp(i,j)+freqsp(k,L) .

compute results(#perm,1) = obsp.
end loop.

* sig levels for the current block of permutations.

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

* two-tailed.
do if (tailed = 2).
compute counter = 0.
loop #i = 1 to nrow(results).
do if   ( sign > 0 and ((results(#i,1) >= obs2) or
                        (results(#i,1) <= obs22))).
compute counter = counter + 1.
else if ( sign < 0 and ((results(#i,1) <= obs2) or
                        (results(#i,1) >= obs22))).
compute counter = counter + 1.
end if.
end loop.
do if (sign ne 0).
compute sigs(#block,1) = counter /  nperms.
end if.
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=(sqrt(cssq (sigs-meansigs)/(nblocks -1)))/(sqrt(nblocks)).
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 meansigs /format="f7.5" /title="Mean Significance Level".
do if (nblocks > 1).
print confid   /title="Percentage for the Confidence Interval:".
print {confidlo,confidhi} /format="f7.5" 
 /title="Low & High Ends of the Confidence Interval".
end if.

end if.
end matrix.