* SEQUENTIAL. * This program computes a variety of basic sequential analysis statistics. 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","Code 16","Code 17","Code 18"}. * Enter the lag number for the analyses here. compute lag = 1. * Can adjacent events be coded the same? Enter "0" if adjacent events can NEVER be the same; Enter "1" if adjacent events can ALWAYS be the same; Enter "2" if some adjacent events can, and others cannot, be the same; then scroll down and enter the appropriate "ONEZERO" matrix for your data. compute adjacent = 1. * Enter "1" for one-tailed tests; enter "2" for two-tailed tests. compute tailed = 2. * 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 }. * transitional frequency matrix. 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. * Warning message for specification error. do if ( adjacent=0 and any(diag(freqs)) = 1). print /title="WARNING". print /title="You have indicated that consecutive or adjacent codes can". print /title="never repeat, (adjacent = 0), yet repeating codes have been". print /title="found in the data. See the main diagonal of the frequency". print /title="matrix. This specification error will result in faulty". print /title="computations for LRX2, z-values, and adjusted residuals.". end if. * initializing. compute lrx2t = {0}. compute rowtots = rsum(freqs). compute coltots = csum(freqs). compute ntrans = msum(rowtots). compute prows = rowtots / ntrans. compute pcols = coltots / ntrans. compute p = make(ncodes,ncodes,-9999). compute et = make(ncodes,ncodes,-9999). compute expfreq = make(ncodes,ncodes,-9999). compute zadjres = make(ncodes,ncodes,-9999). compute pzadjres = make(ncodes,ncodes,1). compute yulesq = make(ncodes,ncodes,-9999). compute var = make(ncodes,ncodes,-9999). compute min = make(ncodes,ncodes,-9999). compute kappa = make(ncodes,ncodes,-9999). compute zkappa = make(ncodes,ncodes,-9999). compute pzkappa = make(ncodes,ncodes,1). compute signs = make(ncodes,ncodes,0). compute n = ntrans + 1. compute nr = rowtots. compute nr(data(nrow(data),1)) = nr(data(nrow(data),1)) + 1. compute pnr = nr / n. loop #i = 1 to ncodes. loop #j = 1 to ncodes. * Note: more refined computations for when adjacent codes cannot repeat appear below, after the above 2 loops are completed. do if (adjacent = 0 and (ntrans - rowtots(#i,1)) > 0 ). compute pcols(1,#j) = coltots(1,#j) / (ntrans - rowtots(#i,1) ). end if. do if (adjacent = 0 and (ntrans - rowtots(#j,1)) > 0 ). compute expfreq(#i,#j) = ( rowtots(#i,1) * coltots(1,#j)) / (ntrans - rowtots(#j,1)). end if. do if (adjacent = 0 and (n - nr(#i,1)) > 0). compute et(#i,#j) = (nr(#i,1) * nr(#j,1)) / (n - nr(#i,1)) . end if. do if (adjacent = 1). compute et(#i,#j) = (nr(#i,1) * nr(#j,1)) / n. compute expfreq(#i,#j) = ( rowtots(#i,1) * coltots(1,#j)) / ntrans . end if. * transitional probabilities. do if (rowtots(#i,1) > 0). compute p(#i,#j) = freqs(#i,#j) / rowtots(#i,1) . end if. * tablewise LRX2. do if ( freqs(#i,#j) > 0 and expfreq(#i,#j) > 0). compute lrx2t=lrx2t+2 * (freqs(#i,#j) * ln(freqs(#i,#j)/expfreq(#i,#j))). end if. * adjusted residuals (z values) & sig levels. do if ( (expfreq(#i,#j)*(1-pcols(1,#j))*(1-prows(#i,1))) > 0). compute zadjres(#i,#j)=(freqs(#i,#j)-expfreq(#i,#j)) /sqrt( expfreq(#i,#j)*(1-pcols(1,#j)) * (1-prows(#i,1)) ). compute pzadjres(#i,#j) = (1 - cdfnorm(abs(zadjres(#i,#j))) ) * tailed. end if. * Yule's Q. compute a = freqs(#i,#j). compute b = rowtots(#i,1) - freqs(#i,#j). compute c = coltots(1,#j) - freqs(#i,#j). compute d = ntrans - rowtots(#i,1) - coltots(1,#j) + freqs(#i,#j). do if ( (a*d + b*c) > 0 ). compute yulesq(#i,#j) = ( a*d - b*c ) / (a*d + b*c). end if. * kappas, z values & sig levels. compute var(#i,#j)=(nr(#i,1)*nr(#j,1)*(n-nr(#j,1))*(n-nr(#i,1)))/ (n**2*(n-1)). do if (var(#i,#j) > 0). compute zkappa(#i,#j) = ( freqs(#i,#j) - et(#i,#j) ) / sqrt(var(#i,#j)). 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. do if ( (min(#i,#j) - et(#i,#j)) ne 0). compute kappa(#i,#j) = ( freqs(#i,#j) - et(#i,#j) ) / ( min(#i,#j) - et(#i,#j) ). do if ( kappa(#i,#j) < 0 ). compute kappa(#i,#j) = ( freqs(#i,#j) - et(#i,#j) ) / et(#i,#j). end if. compute pzkappa(#i,#j) = (1 - cdfnorm(abs(zkappa(#i,#j)))) * tailed. end if. end if. * signs. do if ( freqs(#i,#j) > expfreq(#i,#j) ). compute signs(#i,#j) = 1. else if ( freqs(#i,#j) < expfreq(#i,#j) ). compute signs(#i,#j) = -1. end if. end loop. end loop. do if (adjacent = 0 or adjacent = 2). * maximum likelihood estimation of the expected cell frequencies using iterative proportional fitting (Wickens, 1989, pp. 107-112). compute rsumsf = rsum(freqs). compute csumsf = csum(freqs). compute onezero = make(ncodes,ncodes,1). call setdiag(onezero,0). * the two previous commands create a matrix of ones and zeros that is used in estimating the expected cell frequencies. A "1" indicates that the expected frequency for a given cell is to be estimated, whereas a "0" indicates that the expected frequency for the cell should NOT be estimated, typically because it is a structural zero (codes that cannot follow one another). By default, the matrix that is created by the above commands has zeros along the main diagonal, and ones everywhere else, which will be appropriate for most data sets. However, if your data happen to involve structural zeros that occur in cells other than the cells along the main diagonal, then you must create a ONEZERO matrix with ones and zeros that is appropriate for your data before running any of the commands below. Enter your ONEZERO matrix now. compute expfreq = onezero. loop #ipfloop = 1 to 100. * adjusting by row. compute xr = make(ncodes,1,0). compute rsumse = rsum(expfreq). loop #r = 1 to ncodes. do if ( rsumse(#r,1) > 0). compute xr(#r,1) = rsumsf(#r,1) / rsumse(#r,1). end if. end loop. loop #i = 1 to ncodes. loop #j = 1 to ncodes. do if ( onezero(#i,#j) = 1 ). compute expfreq(#i,#j) = expfreq(#i,#j) * xr(#i,1). end if. end loop. end loop. * adjusting by collumn. compute xc = make(1,ncodes,0). compute csumse = csum(expfreq). loop #c = 1 to ncodes. do if ( csumse(1,#c) > 0). compute xc(1,#c) = csumsf(1,#c) / csumse(1,#c). end if. end loop. loop #i = 1 to ncodes. loop #j = 1 to ncodes. do if ( onezero(#i,#j) = 1 ). compute expfreq(#i,#j) = expfreq(#i,#j) * xc(1,#j). end if. end loop. end loop. compute rdiffs = rsumsf - rsum(expfreq). compute cdiffs = csumsf - csum(expfreq). end loop if ( mmax(rdiffs) < .0001 and mmax(cdiffs) < .0001). print /space = 2. print / title="Maximum likelihood estimation of the expected cell". print / title="frequencies using iterative proportional fitting". do if ( mmax(rdiffs) < .0001 and mmax(cdiffs) < .0001). print #ipfloop/title="converged after the following number of interations:". else. print #ipfloop /title="did NOT converge after the following number of interations:". end if. * tablewise LRX2. compute lrx2t = {0}. loop #i = 1 to ncodes. loop #j = 1 to ncodes. do if ( freqs(#i,#j) > 0 and expfreq(#i,#j) > 0). compute lrx2t=lrx2t+ 2 * (freqs(#i,#j) * ln(freqs(#i,#j)/expfreq(#i,#j))). end if. end loop. end loop. * adjusted residuals for matrices with structural zeros (Christensen, 1997, p. 357). * constructing the design matrix. compute x = make(ncodes**2,1,1). compute y = make(ncodes**2,ncodes-1,0). compute z = make(ncodes**2,ncodes-1,0). loop #i = 1 to ncodes - 1. loop #j = 1 to ncodes. compute y(#i*ncodes+#j,#i) = 1. compute z( (((#j-1)*ncodes)+(#i+1)),#i) = 1. end loop. end loop. compute des1 = { x, y, z }. * pruning values corresponding to cells with structural zeros. compute onezero2 = reshape(onezero,ncodes**2,1). compute dm1 = reshape(expfreq,ncodes**2,1). compute dm2 = make(1,1,-9999). compute des2 = make(1,ncol(des1),-9999). loop #p = 1 to ncodes**2. do if (onezero2(#p,1) = 1). compute dm2 = { dm2 ; dm1 (#p,1) }. compute des2 = { des2 ; des1(#p,:) }. end if. end loop. compute dm2 = dm2(2:nrow(dm2),1). compute des2 = des2(2:nrow(des2),:). compute dm2 = mdiag(dm2). do if ( det(t(des2) * dm2 * des2) ne 0). compute zadjres = make(ncodes,ncodes,0). compute a = des2 * (inv(t(des2) * dm2 * des2)) * t(des2) * dm2. compute acounter = 1. loop #i = 1 to ncodes. loop #j = 1 to ncodes. do if ( onezero(#i,#j) ne 0). compute zadjres(#i,#j) = ( freqs(#i,#j) - expfreq(#i,#j) ) / sqrt( expfreq(#i,#j) * (1 - a(acounter,acounter)) ). compute acounter = acounter + 1. end if. end loop. end loop. else. print /space = 2. print /title="A nonsingular matrix has been identified, which means". print /title="that proper adjusted residuals cannot be computed for". print /title="this data, probably because there are no values for". print /title="one or more codes. Try recoding using sequential". print /title="integers, and redo the analyses. The adjusted". print /title="residuals that are printed below are based on". print /title="equation 5 from Bakeman & Quera (1995, p. 274), and". print /title="are close approximations to the proper values.". print /title="The procedures recommended by Bakeman & Quera". print /title="(1995, p. 276), Haberman (1979), and Christensen". print /title="(1997) cannot be conducted with nonsingular matrices.". end if. loop #i = 1 to ncodes. loop #j = 1 to ncodes. do if ( onezero(#i,#j) = 0). compute zadjres(#i,#j) = 0. compute yulesq (#i,#j) = 0. compute kappa (#i,#j) = 0. compute zkappa (#i,#j) = 0. compute pzadjres(#i,#j) = 1. compute pzkappa (#i,#j) = 1. end if. end loop. end loop. end if. compute b = labels(1,1:ncodes). compute bb = { b,"Totals"}. print /space = 2. print tailed /title="Requested 'tail' (1 or 2) for Significance Tests =". print {freqs, rowtots; coltots, ntrans} /title="Cell Frequencies, Row & Collumn Totals, & N"/cnames=bb/rnames=bb. do if (adjacent = 0 or adjacent = 2). print /space = 1. print /title="The processed ONEZERO matrix appears below.". print /title="In the ONEZERO matrix, a 0 indicates a structural". print /title="zero, and a 1 indicates that an expected cell frequency". print /title="will be estimated.". print onezero /title="ONEZERO matrix:" /cnames=b/rnames=b. end if. print expfreq /format "f5.4" /title="Expected Values/Frequencies" /cnames=b/rnames=b. print p /format "f5.4" /title="Transitional Probabilities" /cnames=b/rnames=b. do if (adjacent = 1). compute df = (ncodes - 1)**2. else. compute df = (ncodes - 1)**2 - (ncodes**2 - msum(onezero)). end if. compute plrx2t = 1 - chicdf(abs(lrx2t),df) . print { lrx2t, df, plrx2t} /format "f12.4" /title="Tablewise Likelihood Ratio (Chi-Square) Test" /cnames={"LRX2","df", "sig."}. print zadjres/format "f6.3" /title="Adjusted Residuals"/cnames=b/rnames=b. print pzadjres /format "f5.4" /title="Significance Levels for the Adjusted Residuals"/cnames=b/rnames=b. print yulesq /format "f6.3" /title="Yule's Q Values"/cnames=b/rnames=b. print kappa /format "f5.4" /title="Unidirectional Kappas"/cnames=b/rnames=b. print zkappa /format "f6.3" /title="z values for the Unidirectional Kappas"/cnames=b/rnames=b. print pzkappa /format "f5.4" /title="Significance Levels for the Unidirectional Kappas" /cnames=b/rnames=b. * Permutation tests of significance. do if (permtest = 1). compute obs2 = reshape(freqs,1,nrow(freqs)*ncol(freqs)). compute obs22 = reshape((expfreq-(freqs-expfreq)),1,nrow(freqs)*ncol(freqs)). compute signs2 = reshape(signs,1,nrow(freqs)*ncol(freqs)). compute sigs = make(nblocks,nrow(freqs)*ncol(freqs),1). 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. compute results(#perm,:) = reshape(freqsp,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 ( results(#i,#j) >= obs2(1,#j) and signs2(1,#j) > 0 ). compute counter = counter + 1. else if ( results(#i,#j) <= obs2(1,#j) and signs2(1,#j) < 0 ). compute counter = counter + 1. 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 ( 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 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 /title="Permutation Tests of Significance:". 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. end matrix.