proc iml; /* Module that returns an "adjustment value" for tied rank. The adjustment value is used in nonparametric tests such as Cuzick's Trend Test.*/ start TiedRankAdj(x); call sortndx(idx, x, 1); sx = x[idx]; /* Adjust for ties (return zero if no ties) */ diff = sx[2:nrow(x)] - sx[1:(nrow(x)-1)]; /* diff consectutive elements */ tieloc = loc(diff=0); tieadj = 0; do while (ncol(tieloc) > 0); tiestart = tieloc[1]; ntied = 1 + sum(sx[tiestart] = sx[tiestart+1:nrow(sx)]); tieadj = tieadj + ntied*(ntied-1)*(ntied+1)/2; tieloc = remove(tieloc, 1:ntied-1); end; return (tieadj); finish; start cuzick(x); /* % CUZICK: Perform the Cuzick's test on trend. % This function provides a Wilcoxon-type test for trend across a group of % three or more independent random samples. % Data must be ordinal and groups are assumed ordered. % % Cuzick J. (1985) "A Wilcoxon-Type Test for Trend," Statistics in Medicine, 4: 87-89. % % This implementation is based on the following: % Cardillo G. (2008) Cuzick’s test: A Wilcoxon-Type Test for Trend % http://www.mathworks.com/matlabcentral/fileexchange/22059 % % SAS/IML implementation by Rick Wicklin, 01DEC2010. % See Cardillo's version for additional details. % See the end of this file for Cardillo's copyright and BSD license. % Syntax: proc iml; d = {0 0 1 1 2 2 4 9 0 0 5 7 8 11 13 23 25 97 2 3 6 9 10 11 11 12 21 0 3 5 6 10 19 56 100 132 2 4 6 6 6 7 18 39 60}; g = j(8,1,1) // j(10,1,2) // j(9,1,3) // j(9,1,4) // j(9,1,5); x = d` || g; run cuzick(x); */ alpha = 0.05; v = x[,1]; /* values */ g = x[,2]; /* groups */ ug = unique(g); k = ncol(ug); /* number of groups */ if k<3 then do; print 'Warning: At least three groups are required'; abort; end; score=1:k; ni = j(1, k); /* allocate vectors */ R = j(1,k); rnk = ranktie(v); /* ranks */ do i = 1 to k; idx = loc(g=ug[i]); ni[i] = ncol(idx); /* number of elements in group */ R[i] = sum(rnk[idx]); /* sum of ranks of each group */ end; T = score * R`; L = score * ni`; N = sum(ni); /* total elements */ ET = L*(N+1)/2; /* mean of T */ labels = {"L" "T" "E(T)"}; print (L || T || ET)[c=labels label="Cuzick's Test For Nonparametric Trend Analysis"]; VarT = (N+1)/12 * (N*score##2*ni` -L##2); /* Variance of T */ z = abs(T-ET)/sqrt(VarT); /* z statistic */ p = 1 - cdf("Normal", z); /* p-value */ labels = {"Var(T)" "z" "1-sided p" "2-sided p"}; print (VarT || z || p || (2*p))[c=labels label="Without Ties Correction"]; if p >= alpha then print 'H_0 accepted.'; else print 'H_0 rejected.'; /* adjust for tied ranks */ AdjVarT = VarT - tiedrankadj(v); z = abs(T-ET)/sqrt(VarT); /* z statistic */ p = 1 - cdf("Normal", z); /* p-value */ print (AdjVarT || z || p || (2*p))[c=labels label="With Ties Correction"]; if p >= alpha then print 'H_0 accepted.'; else print 'H_0 rejected.'; finish; d = {0 0 1 1 2 2 4 9 0 0 5 7 8 11 13 23 25 97 2 3 6 9 10 11 11 12 21 0 3 5 6 10 19 56 100 132 2 4 6 6 6 7 18 39 60}; g = j(8,1,1) // j(10,1,2) // j(9,1,3) // j(9,1,4) // j(9,1,5); x = d` || g; *print x; run cuzick(x); /* Copyright (c) 2008, Giuseppe Cardillo All rights reserved. Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */