/* 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. */ 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; /* Enter/read the data into a matrix with the name "data", as in the following example */ 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,1,1,5,1,6,3,5,3,3,5,3,3,5,3, 6,3,6,6,3,3,6,2,6,3,6,3,3,6,3,3,5,3,6,2,6,2,6,2,6,3,6,6,2,3,6,3,6,2, 6,2,6,2,6,3,3,6,2,6,6,3,3,6,6,3,3,6,2,6,6,4,3,5,3,4,3,6,2,6,1,6,3,6, 3,3,6,1,4,6,3,6,6,3,4,3,6,3,6,3,4,3,5,3,3,6,3,5,3,3,3,6,3,6,3,6,3,6, 3,6,2,6,3,6}; 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; 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; quit;