ࡱ> {}zq` tAbjbjqPqP 8z::t9x1x1x1x12D|C2d2(2222222BBBBBBB$DhGn!C26222626!C226C???26N22B?26B??A/B2X2 52ux1:ABLC0|CB,G>>G/BG/B23?^448222!C!C>"222|C26262626$d  Permutation procd.sas This SAS routine performs a two group permutation test on the Procrustes distance between group mean configurations. In other words, given a distance between the means of two groups, it calculates the probability that both groups derive from the same mean (i.e., that the calculated difference is not significant). Initial group means are first computed from coordinates already aligned by generalized Procrustes analysis. The Procrustes distance is calculated between these means in order to establish the distance benchmark. Next, specimens are randomly assigned to two new groups, new means are calculated, and their Procrustes distance is computed. This randomization is repeated through many iterations (defined by the user) in order to generate a probability distribution of sample mean differences in this combined population. Based on this distribution, one can evaluate the likelihood that the magnitude of difference between the two initial groups could have been sampled from a single population. Note, this routine cannot compute exact permutation tests (something possible with very small group sizes). As presented here, the routine will use permuted group sizes that are identical to the original group sizes. There is an option, however, to use equal sample permuted groups based either on the smaller number of the two original groups, or on a designated group size. Details on this are given below in the annotation for grp1 and grp2 definitions. The section of code commented out (/**/) provides the option of creating a SAS dataset to store all of the permuted Procrustes distances. Removing the bracketing slash-stars will cause this option to execute so these values can be analyzed statistically. Although written for a two-group test, this routine could fairly easily be wrapped in a MACRO (or other) loop in order to compute multiple comparisons. INPUT: One dataset of specimens aligned by GPA, including coordinate data and a single character variable with specimen labels. A character variable specifying group is also necessary in this analysis. In this example, it is generated from the sample data using the specimen labels. However, it can instead be inputted in the same or a separate datafile, requiring only minor modifications to the code. The type of landmark data (2D or 3D) must be specified below (k=); the number of permutation (nreps=) is also designated by the user. OUTPUT: Results are printed to the SAS output window, listing the Procrustes distance between the two groups, the number of permutations used, and the probability that the original Procrustes distance could be sampled from a single combined population. There is also an option to output a SAS dataset containing all of the permuted Procrustes distance. Potential sources of error: Extraneous numeric or character variables left in the dataset when IML begins. Additional numeric variables will usually result in Procrustes distances being too large. Additional character variables will result in an ERROR message at ind=design(SPS) Group variable is numeric instead of character. This will result in an ERROR message at read all var _CHAR_ into SPS. Specimens were not first aligned by GPA. All of the Procrustes distances are computed with an OPA alignment. However, mean configurations need to be computed from GPA aligned specimens. This problem will not be easily spotted. Code annotation: ### indicates lines that need to be changed for different datasets data three; set two; Names dataset three using all of the data from dataset two (here, the sample data) sps=substr(cat,1,2); ### Defines a new variable sps (species) generated by taking the first two characters from the cat variable (i.e., reading a substring of cat starting at character 1 and taking two values). This is the variable which distinguishes the two groups to be tested. drop cat; Deletes the cat variable. As we are testing group differences, the individual specimens labels are not needed. Dropping them from the dataset makes it easier to specify used variable in PROC IML. run; At this point, the dataset needs to be pared down to only two groups of specimens; the only data should be coordinate variables and a single character variable distinguishing the two groups. proc sort data=three; It is important here to sort the data by the group variable to insure that all specimens from the same group are together. by sps; ### Tells SAS which variable to use in sorting. This should be identical to the group variable named in the above data step. run; Executes the sorting process proc means data=three noprint; This procedure will compute the mean configurations for both groups in dataset three. noprint indicates that the means (and many other summary statistics) should not be displayed in the SAS output window. by sps; indicates that separate mean configurations should be computed for each group variable var x1-x60; ### tells SAS the variables for which means should be computed. All coordinate variables should be used. output out=means; output indicates that an output dataset is desired, and out-means designates the name of that dataset means run; Ends the procedure data four; set means; names a new dataset four based on the output from PROC MEANS stored in the means dataset if _STAT_='MEAN'; tells SAS to only include variables in four that are associated with the group means drop _TYPE_ _FREQ_ _STAT_; of those variables, _TYPE_ _FREQ_ _STAT_ are dropped, leaving only the group variable and the mean coordinate variables run; Ends the procedure proc iml; Begins processing commands in Interactive Matrix Language. k=3; ### k is the dimensionality of the data (either 2 or 3). nreps=1000; ### this defines a variable nreps equal to the total number of permutations desired. Change the number accordingly. use three; Opens dataset three for use in IML. This dataset has the aligned coordinates of all specimens. read all into X; Puts all numeric variables into a matrix X. read all var _CHAR_ into SPS; Puts all character variables into SPS. Note, even if the groups used are not species, the name of this matrix should not be changed. use four; read all into MNS; Puts all of the numeric variables from four (i.e., the group mean configurations) into a matrix MNS N=NROW(X); N is the total number of specimens computed as the number of rows in the matrix X L=NCOL(X); L is the number of coordinates, computed as the number of columns in X ind=design(SPS); ind=ind[+,]; ind is the number of specimens in each group. It is calculated using a design matrix of the names stored in SPS and then summing the resulting two rows. p=(L/k); p gives the number of landmarks computed by dividing the number of coordinates (L) by the number of dimensions (k). MN1=MNS[1,]; MN1=shape(MN1,p,k); This section performs an orthogonal Procrustes MN2=MNS[2,]; MN2=shape(MN2,p,k); analysis on the two mean configurations. I=I(p); CaP=j(p,p,1/p); Annotation for this code can be found in the s1=sqrt(trace((I-CaP)*MN1*MN1`*(I-CaP)));  HYPERLINK "http://www.baylor.edu/anthropology/index.php?id=19177" help file for the PROCRUSTES DISTANCE code. s2=sqrt(trace((I-CaP)*MN2*MN2`*(I-CaP))); MN1prime=(I-CaP)*MN1/s1; MN2prime=(I-CaP)*MN2/s2; cov2=MN1prime`*MN2prime; CALL SVD (U,sig,V,cov2); sig2=abs(sig); S=sig/sig2; S=S`; S=diag(S); H=V*S*U`; MN2star=MN2prime*H; y=MN1prime-MN2star; y is a vector of the differences at each coordinate between the superimposed mean configurations d=sqrt(trace(y`*y)); Matrix multiplies the transpose of y by y and then takes the square root of the trace (i.e, the sum of diagonal elements which correspond to the squared differences in each coordinate) of this matrix and defines this distance as d. rho=2*arsin(d/2); Computes Procrustes distance as the angle rho (in radians) as twice the arcsine of d/2. A=SPS[1,]; B=SPS[ind[1,1]+1,]; Retrieves the group names, with A being designated the first group, and B the second. print '...' A 'and ' B ': ' rho; Prints to the SAS output window the Procrustes distance between the group means. count=0; count is a counter variable to index the number of times that Procrustes distance between permuted group means is greater than or equal to the distance between actual group mean (rho) calculated above. small=min(ind); names a variable small which is equal to the number of specimens in the smaller of the two groups. This can be used later to do permutation tests with equal sample sizes. do loop=1 to nreps; sets up a loop for multiple permutations. The variable loop will go from 1 to the number indicated by nreps (defined above), computing new permutations in each iteration. rand=j(N,1,0); defines a new matrix rand using the J function. Here, this creates a matrix of N rows (number of specimens) and 1 column. The default filler for the J function is 1, however, here all cells are filled with 0 rand=normal(rand); redefines rand using the NORMAL function, which generates random numbers. A seed can be placed in the parentheses for random number generation; 0 indicates that the computer clock is used as a seed. When a matrix is named (as in this case), NORMAL will generate a matrix of random numbers equal in rank to that in parentheses. For the seed, it uses whatever value is first in the matrix named in parentheses. As rand was a matrix of zeroes, the clock is used to seed the random number generator. After this step, rand is an N by 1 matrix of random numbers (one for each specimen). ord=rank(rand); creates a matrix ord with all of the ordinal ranks from the random numbers in matrix rand. Note, this does not sort the random numbers, but only indicates the place of each one in numerical sequence. T=j(N,L); defines a matrix T with the J function of N rows (number of specimens) and L columnms (number of coordinates) do i=1 to N; T[i,]=X[ord[i,],]; end; sets up a loop with variable i from one to the number of specimens (N). For each i, the rank from matrix ord is used to indicate the row of the specimen to be added to T. As ord was generated randomly, T is filled in random order with no regard to specimen groups. grp1=T[1:ind[1,1],]; This defines a two new matrices, grp1 and grp2=T[ind[1,1]+1:N,]; grp2, each containing specimens chosen at random (as defined by their order in matrix T). As written here, the number of specimens in grp1 will be equal to the number originally known from the first group, while grp2 will be equal to the number of specimens in the original second group. ### If group sample sizes are unequal but an equal-sample permutation test is desired, replace both instances of ind[1,1] with small. Then replace N with 2*small. This will generate permuted groups equal in size to the smaller of the two group samples. ### Should a specific group size be preferred, replace both instances of ind[1,1] with the number of preferred individuals for each group. Then replace N with twice that number. This option will NOT work properly if the desired group size is more than half of the total sample size. M1=grp1[+,]/NROW(grp1); Computes the mean configurations of the two M2=grp2[+,]/NROW(grp2); permuted groups M1=shape(M1,p,k); This section performs an orthogonal Procrustes M2=shape(M2,p,k); analysis on the two new permuted mean I=I(p); CaP=j(p,p,1/p); configurations. Annotation for this code can be s1=sqrt(trace((I-CaP)*M1*M1`*(I-CaP))); found in the  HYPERLINK "http://www.baylor.edu/anthropology/index.php?id=19177" help file for the PROCRUSTES s2=sqrt(trace((I-CaP)*M2*M2`*(I-CaP))); DISTANCE code. M1prime=(I-CaP)*M1/s1; M2prime=(I-CaP)*M2/s2; cov2=M1prime`*M2prime; CALL SVD (U,sig,V,cov2); sig2=abs(sig); S=sig/sig2; S=S`; S=diag(S); H=V*S*U`; M2star=M2prime*H; y=M1prime-M2star; y is a vector of the differences at each coordinate between the superimposed permuted mean configurations d=sqrt(trace(y`*y)); Matrix multiplies the transpose of y by y and then takes the square root of the trace (i.e, the sum of diagonal elements which correspond to the squared differences in each coordinate) of this matrix and defines this distance as d. rhotemp=2*arsin(d/2); Computes Procrustes distance as the angle rhotemp as twice the arcsine of d/2. if rhotemp >= rho then do; Conditional statement such that every time the Procrustes distance between permuted group means is greater than or equal to the distance between the actual group means, the following statement will be executed. count=count+1; Adds one to the number stored in count. Only occurs if the conditions of the previous statement are met. end; Ends the conditionally executed statements. /* indicates that this section is commented out. It is included here to give the option of examining all of the permuted distances. By removing /* and its paired */ below, a SAS dataset temp will be generated to record all of the permuted mean Procrustes distances. if loop=1 then do; Conditional statement such that the next lines will only be executed if this is the first iteration of the permutations (i.e., loop=1) create temp from rhotemp; creates a SAS dataset temp from the Procrustes distance rhotemp append from rhotemp; adds the value in rhotemp to dataset temp end; Ends the conditionally executed statements if loop>1 then do; Conditional statement such that the next lines will only be executed if this is not the first iteration of the permutations (i.e., loop>1) edit temp; Opens the SAS dataset temp for editing append from rhotemp; Appends the value in rhotemp to the end of the dataset temp end; Ends the conditionally executed statements */ Signifies the end of the section that is commented out end; this ends the loop for each permutation, causing the variable loop to begin again until it has reached the number indicated by nreps alpha=count/nreps; computes the value alpha as the total number of times the Procrustes distance between the two groups was equaled or exceeded by the distances between permuted groups (count), divided by the total number of permutation (nreps) print 'Number of permutations:' nreps; prints in the SAS output window the total number of permutations requested print 'Probability ... same:' alpha; prints in the SAS output window the probability (alpha) that the group means derived from the same population. I.e., the probability that the groups are the same. abort; Ends the IML session.   ~    v w x , 3 c Zqrst '(hiξξʾ־ʾڞ֚š“ h$h$hnh$hS_>*hBch)hFhF>*hFhA>*hAh$hhFhhmuhS_h/6hPCh[h>* h/6>* hBc>*h[hPC>*hK:w x rstlmnopqrs & Fgdngd$ `^``gdS_ `^``gdtAilmnrsĿkP54hB*CJOJQJ^JaJfHphq 4hBcB*CJOJQJ^JaJfHphq :hBc5B* CJOJQJ\^JaJfHphq 4hBcB*CJOJQJ^JaJfHphq hh !h9dhh !hh !hs9!>*hh !h)>* h(*h(*B*CJOJQJ^JaJfHphq 89=DKL!&}~yz{~ 寔vv:hBc5B* CJOJQJ\^JaJfHphq 4h0cB*CJOJQJ^JaJfHphq 4hAB*CJOJQJ^JaJfHphq 4hBcB*CJOJQJ^JaJfHphq 4hBcB*CJOJQJ^JaJfHphq #z{qr  $%jk34P7$8$H$^`Pgd 7$8$H$gdBcP7$8$H$^`Pgd0cpqrv   #%)*-.ijkmnʯʯʯʯvʯvvʯʯ[=:h5B* CJOJQJ\^JaJfHphq 4hB*CJOJQJ^JaJfHphq :hBc5B* CJOJQJ\^JaJfHphq 4hBcB*CJOJQJ^JaJfHphq 4h0cB*CJOJQJ^JaJfHphq 4hBcB*CJOJQJ^JaJfHphq 4hBcB* CJOJQJ^JaJfHphq n4? vvWv<<4h)!B*CJOJQJ^JaJfHphq =h)!h0c>*B*CJOJQJ^JaJfHphq 4h0cB*CJOJQJ^JaJfHphq 4hBcB*CJOJQJ^JaJfHphq :h5B* CJOJQJ\^JaJfHphq 4hB* CJOJQJ^JaJfHphq 4hB*CJOJQJ^JaJfHphq   A F o y ! ! !!l!m!w!!!!!ȭȭȒtV;;;Ȓ4h)!B* CJOJQJ^JaJfHphq :h)!hBcB*CJOJQJ^JaJfHphq :h)!h)!B*CJOJQJ^JaJfHphq 4hBcB*CJOJQJ^JaJfHphq 4h2wB*CJOJQJ^JaJfHphq 4h)!B*CJOJQJ^JaJfHphq 7h)!>*B*CJOJQJ^JaJfHphq  ! !l!m!!!""##U###$$$$$%1%B%L%a%b% 7$8$H$gdBcP7$8$H$^`Pgd)!!""# # #&#T#]#^#u######$$$W$X$Y$ʬʑʬʑʬʑʑrMrIjh2wh&.B*CJOJQJU^JaJfHphq =jh2wB*CJOJQJU^JaJfHphq 4h2wB*CJOJQJ^JaJfHphq :hBc5B* CJOJQJ\^JaJfHphq 4hBcB*CJOJQJ^JaJfHphq 4h)!B*CJOJQJ^JaJfHphq Y$b$c$$$$$%u%%%%%%&&&&&&&ŪtY;;:hBc5B* CJOJQJ\^JaJfHphq 4h2wB* CJOJQJ^JaJfHphq 4hBcB*CJOJQJ^JaJfHphq 4hBcB*CJOJQJ^JaJfHphq 4h2wB*CJOJQJ^JaJfHphq =jh2wB*CJOJQJU^JaJfHphq 5h2wh2w0JCJOJQJ^JaJfHq b%%%&&O'P'''@(A(B())))**++--.P7$8$H$^`Pgdl4P7$8$H$^`Pgd{a 7$8$H$gdBcP7$8$H$^`Pgd2w&'O'P'V'W'e'f'g'h'j'k'o''''''''''''?(ʯvvvv[@@@4hBcB* CJOJQJ^JaJfHphq 4h{aB* CJOJQJ^JaJfHphq :hBc5B* CJOJQJ\^JaJfHphq 4hBcB*CJOJQJ^JaJfHphq 4hBcB* CJOJQJ^JaJfHphq 4h2wB*CJOJQJ^JaJfHphq 4h{aB*CJOJQJ^JaJfHphq ?(H(I(J()')))))******d++++--..>/?/@/E/F/`/c/d//y000嬑v[4h`TJB*CJOJQJ^JaJfHphq 4hBcB*CJOJQJ^JaJfHphq 4hl4B*CJOJQJ^JaJfHphq 4h{aB*CJOJQJ^JaJfHphq :hBc5B* CJOJQJ\^JaJfHphq 4hBcB*CJOJQJ^JaJfHphq #..?/@/z0{0013"4#4j44445b5536J6a6x666 7$8$H$^gdr<P7$8$H$^`Pgd`TJP7$8$H$^`Pgdl4 7$8$H$gdBc000000000000000011\2`2o2y233:3嬑嬑u\A4hFB*CJOJQJ^JaJfHphq 1h`TJh`TJCJOJQJ^JaJfHq 7h`TJ>*B*CJOJQJ^JaJfHphq 4hr<B*CJOJQJ^JaJfHphq 4h`TJB*CJOJQJ^JaJfHphq :hBc5B* CJOJQJ\^JaJfHphq 4hBcB*CJOJQJ^JaJfHphq :3>3M3W33!4"4:4i44444445)5ȯȔy^y^y^Cy^y4hr<B* CJOJQJ^JaJfHphq 4hr<B*CJOJQJ^JaJfHphq 4hBcB*CJOJQJ^JaJfHphq 4hB*CJOJQJ^JaJfHphq 1h`TJhFCJOJQJ^JaJfHq 4hFB*CJOJQJ^JaJfHphq 7hF>*B*CJOJQJ^JaJfHphq )5*5.5a5555555555#626x6|66X7Y7ǬǬhMǬ2Ǭ4hBcB*CJOJQJ^JaJfHphq 5h2whr<0JCJOJQJ^JaJfHq Ijh2whr<B*CJOJQJU^JaJfHphq =jhr<B*CJOJQJU^JaJfHphq 4hr<B*CJOJQJ^JaJfHphq 4hBcB*CJOJQJ^JaJfHphq :hBc5B* CJOJQJ\^JaJfHphq 66666Y7Z7]8^888990:1:h:i:j:{;|;<<y<z<<P7$8$H$^`Pgd4P7$8$H$^`Pgdr< 7$8$H$gdBcY7\8]8^8f8g8p8q8s88899999/:1:4:5:\:g:i:j:l:ʬʬʬv[@4hBcB* CJOJQJ^JaJfHphq 4h4B* CJOJQJ^JaJfHphq 4hBcB*CJOJQJ^JaJfHphq 4h4B*CJOJQJ^JaJfHphq :hBc5B* CJOJQJ\^JaJfHphq 4hBcB*CJOJQJ^JaJfHphq 4hr<B*CJOJQJ^JaJfHphq l:m:n:; ;_;d;y;z;{;|;;<<<4<x<y<z<<<<<<<<<==T=W=X=[==========>϶϶϶϶϶l϶϶϶1h4hBcCJOJQJ^JaJfHq 4hBcB* CJOJQJ^JaJfHphq +h4CJOJQJ^JaJfHq 1h4h4CJOJQJ^JaJfHq +h4CJOJQJ^JaJfHq 4h4B* CJOJQJ^JaJfHphq *<<<<====>>O>P>>>??@@@@QARAtA 7$8$H$gd(>>#>N>O>P>R>S>>>>>>?,?k?˰˰˚fKf0f04h4B*CJOJQJ^JaJfHphq 4hBcB*CJOJQJ^JaJfHphq 4hBcB*CJOJQJ^JaJfHphq 1h4hBcCJOJQJ^JaJfHq +h4CJOJQJ^JaJfHq 4hBcB* CJOJQJ^JaJfHphq 4h4B* CJOJQJ^JaJfHphq 1h4h4CJOJQJ^JaJfHq k?|??@@@5@<@@@@@@@@@PAXAsAtAʯ寔y[:h(*B*pht9zwxrst   l m n o p q r s 67WX)*  -.~z{qr  $%jk34  lmU1BLabOP@ A B !!!!""##%%&&?'@'z({(()+",#,j,,,,-b--3.J.a.x.......Y/Z/]0^000110212h2i2j2{3|344y4z44444555566O6P666778888Q9R9v900000000000000 0 0 000000000000000000000000000000000000000000000@0@000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000wxrt   l m r 67WX)*  -.~z{qr  $%jk34  lmU1BLabOP@ A B !!!!""##%%&&?'@'z({(()+",#,j,,,,-b--3.J.a.x.......Y/Z/]0^000110212h2i2j2{3|344y4z44444555566O6P666778888Q9R9v90000000000000 0 0 00000000000000000000000000000000000000@0@0@0@000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000i8n !Y$&?(0:3)5Y7l:>k?tA!$&'(*+,./1245689;<sb%.6<tA"%)-037:tA#Xb---t9XX5ڒt!6ڒ<#7ڒ$#8ڒ\Q9ڒN:ڒ";ڒ^Q<ڒ$!=ڒToQ>ڒR?ڒQ@ڒ܄QAڒBڒgQCڒT MDڒl MEڒ|#Fڒ/RGڒdQ aa$$((((o*o*M+M+v9       dd$$((((r*r*P+P+v9 8 *urn:schemas-microsoft-com:office:smarttagsCity9*urn:schemas-microsoft-com:office:smarttagsState9*urn:schemas-microsoft-com:office:smarttagsplace L  g l  r u s w x 7:;AilDK"%*- ]^cdghk"9=OR #ad!!"!%!!!S"X"%%&&'&'C'D'O'P'U'X'Y'Z'''''''((((o*r*M+P+-"-e-i-s-v----. ....>.A.U.X.....\/`///^0e0h0m0000000+424p4w44444555577&7+7 8868;8v9 X\*- .2!{~ rv %)47 Y]RV!!""""T#U###%%&&@'B'((((o*s*M+Q+&,+,m,r,,,,,--e-j--.:.<.Q.S.....\/a/^0f0111242|3~34 4z44444455556"6S6\688889P9R9W9v9333333333333333333333333333333333333333333333333333333333333333333333333333333rs   ' (  k m q r s #cs9s9v9s w cv9P8u^`OJPJQJ^Jo(-^`OJQJ^Jo(hHopp^p`OJQJo(hH@ @ ^@ `OJQJo(hH^`OJQJ^Jo(hHo^`OJQJo(hH^`OJQJo(hH^`OJQJ^Jo(hHoPP^P`OJQJo(hHP8u$:        -,[0$ 0ch !)!s9!*0l4/6r<E?PC(1TableFGWordDocument8zSummaryInformation(jDocumentSummaryInformation8rCompObjq  FMicrosoft Office Word Document MSWordDocWord.Document.89q