/*  TWOCELLS - 2 */

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



options nocenter nodate nonumber linesize=90;   title;

proc iml;

/* Enter the number of possible code values here. */
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. */
i =  1;
j =  6;
k =  5;
L =  2;

/* Enter labels for the codes (5 characters maximum), if desired, here. */
labels={Code1 Code2 Code3 Code4   Code5  Code6  Code7
        Code8 Code9 Code10 Code11 Code12 Code13 Code14 Code15};

/* Enter the lag number for the analyses here. */
lag = 1;

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

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

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

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

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

/* Enter "95" for 95% confidence intervals;
   enter "99" for 99% confidence intervals */
confid = 95;

/* Set 'seed' to a non-negative integer; 1953125 is good. */
seed = 1953125;




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


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 };



dyadnum   = alldata[1,1];
data = { 0 };
out = j(1,5,0);

/* setting up the data for each group/dyad. */
do dydloop = 1 to nrow(alldata);
if ( dydloop > 1 & alldata[dydloop,1] < 1000 ) then;
data = ( data // alldata[dydloop,1] );

if ( dydloop > 1 & alldata[dydloop,1] >= 1000 )  then do;
data = data[2:nrow(data), 1];



/* start of group/dyad analyses. */





freqs = j(ncodes,ncodes,0);
do c = 1 to nrow(data);
if ( c + lag <= nrow(data) )    then
freqs[(data[c,1]),(data[(c+lag),1])] =
        freqs[(data[c,1]),(data[(c+lag),1])]  + 1;
end;

rowtots =j(nrow(freqs),1);
coltots =j(1,ncol(freqs));
do t = 1 to ncodes;
dumr = freqs[t,];
rowtots[t,1]=sum(dumr);
dumc = freqs[,t];
coltots[1,t]=sum(dumc);
end;
n = nrow(data);
nr  = rowtots;
nr[data[n,1]] = nr[data[n,1]] + 1;
ett = (nr[i,1]*nr[j,1] + nr[k,1]*nr[L,1] ) /n;
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));
zkappa  = ( (freqs[i,j]+freqs[k,L]) - ett ) / sqrt(var);
pzkappa = (1 - probnorm(abs(zkappa))) * tailed;
if (nr[i,1] <= nr[j,1])    then  minij = nr[i,1];
else                             minij = nr[j,1];
if (nr[k,1] <= nr[L,1])    then  minkL  = nr[k,1];
else                             minkL  = nr[L,1];
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));
if ( kappa < 0)    then
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));


print, "Requested 'tail' (1 or 2) for Significance Tests =", tailed;
print, "Session/Dyad/Group Number:", dyadnum;
b = labels[1,1:ncodes];
bb = ( b || {Totals});
transf = ( (freqs || rowtots) // (coltots || sum(rowtots)) );
print, "Cell Frequencies, Row & Collumn Totals, & N" ,
         transf [rowname=bb colname=bb format=7.0];
ijkl = ( i || j || k || L );
print, "Simultaneous Two-Cell Test for the following values of i, j, k, & L",
        ijkl[colname={I J K L } format=12.4];
stats = (ett || (freqs[i,j]+freqs[k,L]) || kappa || zkappa || pzkappa );
print, "Results:", stats[colname={Expected Observed kappa z sig} format=12.4];




/* Permutation tests of significance */
if (permtest = 1) then do;

obs2    = freqs[i,j]+freqs[k,L];
obs22   = ett-(obs2-ett);
sign = 0;
if (kappa > 0)  then;
sign =  1;
if (kappa < 0)  then;
sign = -1;
else;
sign = 0;
sigs = j(nblocks,1,1);

do block = 1 to nblocks;
print, "Currently computing for block :", block;

results = j(nperms,1,-9999);

do perm = 1 to nperms;

/* permuting the sequences; alogrithm from Castellan 1992. */

/* when adjacent codes may be the same */
datap = data;
if (adjacent = 1)  then do;
do ii = 1 to (nrow(datap) -1);
kay = int( (nrow(datap) - ii + 1) * uniform(seed) + 1 )  + ii - 1;
d = datap[ii,1];
datap[ii,1] = datap[kay,1];
datap[kay,1] = d;
end;
end;

/* when adjacent codes may NOT be the same */
if (adjacent = 0) then do;
datap = ( 0 // data // 0 );
do ii = 2 to (nrow(datap) - 2);
limit = 10000;
do jj = 1 to limit;
kay = int( ((nrow(datap)-1) - ii + 1) * uniform(seed) + 1 )  + ii - 1;
if ( (datap[ii-1,1]  ^= datap[kay,1]) & (datap[ii+1,1]  ^= datap[kay,1])
   & (datap[kay-1,1] ^= datap[ii,1])  & (datap[kay+1,1] ^= datap[ii,1]) ) then;
jj = limit;
end;
d = datap[ii,1];
datap[ii,1] = datap[kay,1];
datap[kay,1] = d;
end;
datap = datap[2:(nrow(datap)-1),];
end;

/* transitional frequency matrix for permuted data. */
freqsp = j(ncodes,ncodes,0);
do c = 1 to nrow(datap);
if ( c + lag <= nrow(datap) )    then
freqsp[(datap[c,1]),(datap[(c+lag),1])] =
 freqsp[(datap[c,1]),(datap[(c+lag),1])]  + 1;
end;

/* two-cell frequency for permuted data */
obsp = freqsp[i,j]+freqsp[k,L];

results[perm,1] = obsp;
end;


/* sig levels for the current block of permutations */

/* one-tailed */
if (tailed = 1) then do;
counter = 0;
do ii = 1 to nrow(results);
if   ( results[ii,1] >= obs2 & sign > 0 )    then;
counter = counter + 1;
if   ( results[ii,1] <= obs2 & sign < 0 )    then;
counter = counter + 1;
end;
if (sign ^= 0) then;
sigs[block,1] = counter / nperms;
end;

/* two-tailed */
if (tailed = 2) then do;
counter = 0;
do ii = 1 to nrow(results);
if ( sign > 0 &  ((results[ii,1] >= obs2) |
                  (results[ii,1] <= obs22))) then;
counter = counter + 1;
if ( sign < 0 &  ((results[ii,1] <= obs2) |
                  (results[ii,1] >= obs22))) then;
counter = counter + 1;
end;
if (sign ^= 0) then;
sigs[block,1] = counter /  nperms;
end;

end;


/* mean significance levels and confidence intervals */
if  (confid = 95 & tailed = 1)  then
z = 1.645;
if  (confid = 95 & tailed = 2)  then
z = 1.96;
if  (confid = 99 & tailed = 1)  then
z = 2.326;
if  (confid = 99 & tailed = 2)  then
z = 2.576;

meansigs = sum(sigs) / nblocks;
if (nblocks > 1) then do;
semeans=(sqrt(ssq(sigs-meansigs)/(nblocks-1)))/(sqrt(nblocks));
confidhi = meansigs + z # semeans;
confidlo = meansigs - z # semeans;
end;

print, "Number of permutations per block: ", nperms;
print, "Number of blocks of permutations: ", nblocks;
print, "Mean Significance Level", meansigs[format=7.5];
if (nblocks > 1) then do;
print, "Percentage for the Confidence Intervals:", confid;
print, "High End of the Confidence Interval", confidhi;
print, "Low End of the Confidence Interval", confidlo;
end;

end;


/* Would you like to save any matrix output data for use in SAS
  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 CREATE command.  You must first provide 
  a valid SAS-data-set name for your system.  The command is currently 
  set to save the kappa output matrix.  Substitute an alternative
  matrix name of your choice. */

matrxout = kappa;


out=(out // (dyadnum || ett || (freqs[i,j]+freqs[k,L]) || kappa || zkappa));
data    = {0};
dyadnum   = alldata[dydloop,1];
end;
end;
out = out[2:nrow(out),];


/* The following CREATE command will save the above-specified matrix
   output data in a SAS file that can be accessed for further processing.
   You must provide a legal file name for your system and then "activate"
   the CREATE command, which is currently  a mere comment. */

create dataset from out;
append from out;

quit;


proc means data=dataset;
run;