/* NONPARADOM */ /* This program tests for nonparallel dominance, a form of asymmetry in predictability, between i to j, and k to L, 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. */ options nocenter nodate nonumber linesize=90; title; proc iml; /* Enter the number of possible code values here. */ ncodes = 6; /* This program tests the difference in predictability between i to j, and k to L Enter the code values for your desired test below. */ i = 5; j = 3; k = 6; 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; pd = j(ncodes,ncodes,-9999); et = j(ncodes,ncodes,-9999); 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; prow = nr / (sum(nr)); do ii = 1 to ncodes; do jj = 1 to ncodes; if (nr[ii,1] > 0 & nr[jj,1] > 0) then pd[ii,jj] = ( (freqs[ii,jj] / nr[ii,1]) - prow[jj,1] ) / prow[jj,1]; if (nr[ii,1] > 0 & nr[jj,1] > 0) then et[ii,jj] = (nr[ii,1] * nr[jj,1]) / n; end; end; 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]; /* kappas. */ if ( freqs[i,j] = et[i,j] ) then do; kappa = 0; case = 0; end; /* Wampold's 1st case. */ if ( freqs[i,j] > et[i,j] & freqs[k,L] >= et[k,L] ) then do; kappa=(nr[k,1]*nr[L,1]*freqs[i,j]-nr[i,1]*nr[j,1]*freqs[k,L])/ ( nr[k,1]*nr[L,1]*minij - (nr[i,1]*nr[j,1]*nr[k,1]*nr[L,1]/n)); if ( kappa < 0 ) then kappa=(nr[k,1]*nr[L,1]*freqs[i,j]-nr[i,1]*nr[j,1]*freqs[k,L])/ ( nr[i,1]*nr[j,1]*minkL - (nr[i,1]*nr[j,1]*nr[k,1]*nr[L,1]/n)); case = 1; end; /* Wampold's 2nd case. */ if ( freqs[i,j] < et[i,j] & freqs[k,L] <= et[k,L] ) then do; kappa=(nr[k,1]*nr[L,1]*freqs[i,j]-nr[i,1]*nr[j,1]*freqs[k,L])/ (-1*(nr[i,1]*nr[j,1]*nr[k,1]*nr[L,1]/n)); case = 2; end; /* Wampold's 3rd case. */ if ( freqs[i,j] > et[i,j] & freqs[k,L] <= et[k,L] ) then do; kappa=(nr[k,1]*nr[L,1]*freqs[i,j]+nr[i,1]*nr[j,1]*freqs[k,L]- 2*(nr[i,1]*nr[j,1]*nr[k,1]*nr[L,1]/n) )/ ( nr[k,1]*nr[L,1]*minij - (nr[i,1]*nr[j,1]*nr[k,1]*nr[L,1]/n)); if ( kappa < 0 ) then kappa=(nr[k,1]*nr[L,1]*freqs[i,j]+nr[i,1]*nr[j,1]*freqs[k,L]- 2*(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); case = 3; end; /* Wampold's 4rth case. */ if ( freqs[i,j] < et[i,j] & freqs[k,L] >= et[k,L] ) then do; kappa=(nr[k,1]*nr[L,1]*freqs[i,j]+nr[i,1]*nr[j,1]*freqs[k,L]- 2*(nr[i,1]*nr[j,1]*nr[k,1]*nr[L,1]/n) )/(- 1*(nr[i,1]*nr[j,1]*nr[k,1]*nr[L,1]/n)); if (kappa < 0) then kappa=(nr[k,1]*nr[L,1]*freqs[i,j]+nr[i,1]*nr[j,1]*freqs[k,L]- 2*(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) - nr[i,1]*nr[j,1]*minkL ); case = 4; end; /* observed frequency, expected frequency, variance, & z. */ zeqk = 9999; /* same direction. */ if ( (pd[i,j]>=0 & pd[k,L]>=0) | (pd[i,j]<=0 & pd[k,L]<=0) ) then do; obs=nr[k,1]*nr[L,1]*freqs[i,j]-nr[i,1]*nr[j,1]*freqs[k,L]; ett = 0; var=(nr[i,1]*nr[j,1]*nr[k,1]*nr[L,1]* (n*nr[i,1]*nr[j,1] + n*nr[k,1]*nr[L,1] - nr[i,1]*nr[j,1]*nr[k,1] - nr[i,1]*nr[j,1]*nr[L,1] - nr[i,1]*nr[k,1]*nr[L,1] - nr[j,1]*nr[k,1]*nr[L,1] ))/ (n*(n-1)); zkappa = obs / sqrt(var); zeqk = zkappa; end; /* different directions. */ if ( (pd[i,j]<=0 & pd[k,L]>=0) | (pd[i,j]>=0 & pd[k,L]<=0) ) then do; obs=nr[k,1]*nr[L,1]*freqs[i,j]+nr[i,1]*nr[j,1]*freqs[k,L]; ett = 2*nr[i,1]*nr[j,1]*nr[k,1]*nr[L,1]/n; var=(nr[i,1]*nr[j,1]*nr[k,1]*nr[L,1]* (4*nr[i,1]*nr[j,1]*nr[k,1]*nr[L,1] + n##2 *nr[i,1]*nr[j,1] + n##2 *nr[k,1]*nr[L,1] - n*nr[i,1]*nr[j,1]*nr[k,1] - n*nr[i,1]*nr[j,1]*nr[L,1]- n*nr[i,1]*nr[k,1]*nr[L,1] - n*nr[j,1]*nr[k,1]*nr[L,1] )) / (n##2 *(n-1)); zkappa = ( obs - ett ) / sqrt(var); if ( ((pd[i,j] = 0 | pd[j,i] = 0)) & zeqk < zkappa ) then zk = zeqk; end; pzkappa = (1 - probnorm(abs(zkappa))) * tailed; sign = 0; if ( kappa > 0 ) then sign = 1; if ( kappa < 0 ) then sign = -1; 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]; print, "Expected Values/Frequencies", et[rowname=b colname=b format=5.4]; ijkl = ( i || j || k || L ); print, "Nonparallel Dominance Test for the following values of i, j, k, & L", ijkl[colname={I J K L= } format=12.4]; print, "Nonparallel Dominance 'Case' Type (Wampold, 1989):", case; print, "Case 1: i increases j, and j increases i"; print, "Case 2: i decreases j, and j decreases i"; print, "Case 3: i increases j, and j decreases i"; print, "Case 4: i decreases j, and j increases i"; stats = ( ett || obs || kappa || (zkappa#sign) || pzkappa ); print, "Results:", stats[colname={expected observed kappa z sig} format=12.4]; print, "The above Expected & Observed values are weighted; see Wampold, 1989, p 184"; /* Permutation tests of significance */ if (permtest = 1) then do; obs2 = obs; obs22= ett-(obs2-ett); 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; /* nonparallel dominance frequency for permuted data */ np = nrow(datap); nrp = j(ncodes,1,0); do t = 1 to ncodes; dumr = freqsp[t,]; nrp[t,1]=sum(dumr); end; nrp[datap[np,1]] = nrp[datap[np,1]] + 1; if (case = 1 | case = 2) then; obsp=nrp[k,1]*nrp[L,1]*freqsp[i,j]-nrp[i,1]*nrp[j,1]*freqsp[k,L]; if (case = 3 | case = 4) then; obsp=nrp[k,1]*nrp[L,1]*freqsp[i,j]+nrp[i,1]*nrp[j,1]*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 ( case = 1 | case = 3 ) then do; if ( sign > 0 & results[ii,1] >= obs2 ) then; counter = counter + 1; if ( sign < 0 & results[ii,1] <= obs2 ) then; counter = counter + 1; end; if ( case = 2 | case = 4 ) then do; if ( sign > 0 & results[ii,1] <= obs2 ) then; counter = counter + 1; if ( sign < 0 & results[ii,1] >= obs2 ) then; counter = counter + 1; end; 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 ( case = 1 | case = 3 ) then do; 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 ( case = 2 | case = 4 ) then do; 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; 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;