sil=all ; screen 39x80 ; if ( ntrees < 0 ) errmsg Must read data (and trees) before running ; end macro *10 ( ( ( 2*ntax+2 ) * 180 ) + 1000 ) ; macro [ ( ( ( ( ( 2*ntax+2 ) * 400 ) * 4 + 1000 ) / 1000 ) + 5000 ) ; macro=; /*** This script is a refinement of the original script of Giannini & Goloboff (2010), delcor.run. It is intended for use in large data sets where there are too many reconstructions to be handled by the original script. This version also adds the option of saving the reconstructions for later use (so that they do not have to be generated again), and the option of using fewer significant digits and/or averages for the terminal taxa (which may speed up enumeration of reconstructions). The original script is retained for documentation purposes. ***/ macfloat 6 ; /************************************************************\ /************************************************************\ | | | M A I N B O D Y O F S C R I P T | | | | This contains the main set of instructions, and | | (with "goto") calls the rest of the routines | | that are used in the script. | | | \************************************************************/ var: numsigdigs averagit rdsvtype savefilename[ 100 ] signif matches logged dograph userlengths reccycles numstordx numstordy notest reclo rechi testdown maxstorrays twotailed incx [ (2*ntax+2) ] incy [ (2*ntax+2) ] visitd [ (2*ntax+2) ] distosrc [ (2*ntax+2) ] radius randwts lor hir reconum sample numrepls outnod innod tmpradius usertree ispos userx worsr usery minincx minincy myx myy drag [ (2*ntax+2) ] xpairs [ (2*ntax+2) ] ypairs [ (2*ntax+2) ] blength [ (2*ntax+2) ] numpairs wts [ (2*ntax+2) ] delfac didpair stasis dax day i ; /******************** Explanation of main variables used in script: dograph : when set to N, plots pairs for set of reconstructions number N maxstorrays : maximum number of reconstructions to store for X and Y stordx[] & stordy[] : stored arrays for increments (to facilitate randomization, and to make sure randomization uses a random sample of reconstructions) notest : skip the randomization testdown : test as tree-down. That is, check whether change Y is always preceded by a change in X. Alternative (default) is checking whether a change in X always implies a change in Y. twotailed : make test two-tailed radius : how many nodes to go up (or down) to search for a match randwts : whether to use a randomized radius value reccycles : if using randwts, how many cycles to do for observed correlation lor & hir : value of lowest/highest correlation reclo & rechi : number of reconstructions which produces lowest/highest correlation usertree : number of tree to use for testing userx & usery : number of characters, X and Y, to test minincx & minincy : threshold of change in X and Y to be considered blength [] : array to store branch lengths provided by user userlengths : use branch lengths provided by user (otherwise, all unity) delfac : delay factor, in formula to set weight as 1 / ( 1 + ( Dist x delfac ) ) stasis : **********/ goto = %0 ; goto READPARMS %1 %2 %3 %4 %5 %6 %7 %8 %9 %(10) %(11) %(12) %(13) %(14) %(15) %(16) %(17) %(18) %(19) %(20) ; goto CHECKPARMS ; set maxstorrays 'sample' ; var: stordx [ 'maxstorrays' (2*ntax+2) ] stordy [ 'maxstorrays' (2*ntax+2) ] ; report- ; lquote=; lquote[ ; naked=; if ( 'userlengths' ) quote &10SETTING USER BRANCH LENGTHS ...&10 ; /********************************************************** Read user-supplied branch lengths, when so requested */ loop 0 nnodes['usertree'] if ( #1 == root ) continue ; end if ( eqstring[ x$ttag #1 x ] ) errmsg Sorry, branch nr. #1 has no length defined... &10; end set blength[#1] ( $ttag #1 ) ; stop quote &32 ...all branch lengths set successfully.&10; else /************************************************************* When no user-supplied branch lengths, all branches unity */ loop 0 nnodes['usertree'] set blength[#1] 1 ; stop end if ( !exstatus ) proc/; end if ( !'randwts' ) set reccycles 1 ; end if ( !'testdown' ) goto CORRELATE_UP ; else goto CORRELATE_DOWN ; end if ( 'dograph' || 'notest' ) proc/; end if ( 'numrepls' == 0 ) proc/; end if ( ( 'lor' <= 0 ) && ( 'hir' >= 0 ) ) sil - all ; quote &10NO NEED TO RANDOMIZE ANYTHING!!&10 ; proc/; end if ( 'lor' > 0 ) set ispos 1 ; else set ispos 0 ; end if ( !'testdown' ) goto RANDOMIZE_UP ; else goto RANDOMIZE_DOWN ; end proc/; /************************************************************\ /************************************************************\ | | | I N I T I A L I Z E A N D | | R E A D P A R A M E T E R S | | | \************************************************************/ label READPARMS ; set usertree 0 ; /* the tree for the test */ set userx 0 ; /* primary (X) variable */ set usery 1 ; /* secondary (Y) variable */ set minincx 0.0001 ; /* minimum change in X */ set minincy 0.0001 ; /* minimum change in Y */ set radius 4 ; /* max. nr. of nodes to visit seeking a change in y */ set numpairs (-1) ; /* initial number of XY pairs */ set delfac 1 ; /* factor to downweight pairs, wt = 1 / ( 1 + ( B*factor) ), where B=branches in between */ set stasis 0 ; /* are cases of (0,0) or (0,1) considered? */ set randwts 0 ; /* randomize weights for cases of stasis? */ set sample 50 ; /* number of combinations of reconsts. to look at */ set numrepls 100 ; /* num replications */ set dograph 0 ; /* graphically display results and quit */ set testdown 0 ; /* does Y change only if X changed? */ set reccycles 5 ; /* number of cycles for each reconst (randomized radius) */ set userlengths 0 ; /* Dont use supplied branch lengths */ set twotailed 0 ; /* Make test two-tailed */ set notest 0 ; /* notest: skip randomization */ set numsigdigs 3 ; /* number of significant digits to use */ set averagit 0 ; /* use average of min/max for terminals */ set rdsvtype 0 ; /* dont use files; 1 save to a file; 2 read from file */ var: i j k ; if ( argnumber == 1 ) if ( eqstring [ %1 help ] ) goto DISPLAY_HELP ; return 0 ; end end if ( argnumber ) loop 1 argnumber set i #1 ; set j #1 + 1 ; set k #1 + 2 ; if ( eqstring [ %('i') tree ] ) set usertree %('j') ; var usertree ; setloop 'k' ; else if ( eqstring [ %('i') chars ] ) set userx %('j') ; var userx ; set usery %('k') ; var usery ; setloop ( #1 + 3 ) ; else if ( eqstring [ %('i') minincx ] ) set minincx %('j') ; var minincx ; setloop ( #1 + 2 ) ; else if ( eqstring [ %('i') minincy ] ) set minincy %('j') ; var minincy ; setloop ( #1 + 2 ) ; else if ( eqstring [ %('i') radius ] ) set radius %('j') ; var radius ; setloop ( #1 + 2 ) ; else if ( eqstring [ %('i') stasis ] ) set stasis 1 ; var stasis ; else if ( eqstring [ %('i') delfac ] ) set delfac %('j') ; var delfac ; setloop ( #1 + 2 ) ; else if ( eqstring [ %('i') sample ] ) set sample %('j') ; var sample ; setloop ( #1 + 2 ) ; else if ( eqstring [ %('i') cycles ] ) set reccycles %('j') ; var reccycles ; setloop ( #1 + 2 ) ; else if ( eqstring [ %('i') repls ] ) set numrepls %('j') ; var numrepls ; setloop ( #1 + 2 ) ; else if ( eqstring [ %('i') dograph ] ) set dograph %('j') ; var dograph ; setloop ( #1 + 2 ) ; else if ( eqstring [ %('i') twotailed ] ) set twotailed 1; var twotailed ; else if ( eqstring [ %('i') userlengths ] ) set userlengths 1; var userlengths ; else if ( eqstring [ %('i') testdown ] ) set testdown 1; var testdown ; else if ( eqstring [ %('i') notest ] ) set notest 1; var notest ; else if ( eqstring [ %('i') randwts ] ) set randwts 1; var randwts ; else if ( eqstring [ %('i') saveto ] ) set rdsvtype 1 ; set savefilename $ %('j'); setloop ( #1 + 2 ) ; else if ( eqstring [ %('i') rdfrom ] ) set rdsvtype 2 ; set savefilename $ %('j'); setloop ( #1 + 2 ) ; else if ( eqstring [ %('i') average ] ) set averagit 1 ; var averagit ; else if (eqstring [ %('i') sigdigs ] ) set numsigdigs %('j') ; var numsigdigs ; setloop ( #1 + 2 ) ; else errmsg Wrong argument: %('i'); end end end end end end end end end end end end end end end end end end end end stop else if ( !windows ) return 1 ; end var : stgetminx[4] stgetminy[4] stgetradius[4] stgetdelfac[4] stgetnumrepls[4] havtreetags ; set stgetminx $'/.3minincx' ; set stgetminy $'/.3minincy' ; set stgetradius $'/.3radius' ; set stgetdelfac $'/.3delfac' ; set stgetnumrepls $'/.0numrepls' ; set havtreetags ( eqstring[ x$ttag 0 x ] == 0 ) ; opendlg 15 15 570 360 Measure and test correlation... frame 10 10 160 120 Choose chars. and tree ; spin 0 ntrees usertree 15 30 75 20 tree nr.; spin 0 nchar userx 15 55 75 20 char. X; spin 0 nchar usery 15 80 75 20 char. Y; subdlg 50 105 110 20 H e l p ; 10 10 220 140 Choose chars. and tree ; showtxt 10 10 190 60 @@ CHOOSE_CHARS ; closedlg frame 10 140 160 145 General settings ; gettxt stgetminx 15 160 50 20 ; showtxt 75 160 100 18 Min.incr.X; gettxt stgetminy 15 185 50 20 ; showtxt 75 185 100 18 Min.incr.Y; gettxt stgetradius 15 210 50 20 ; showtxt 75 210 100 18 (Max)radius; gettxt stgetdelfac 15 235 50 20 ; showtxt 75 235 100 18 Delay factor; subdlg 50 260 110 20 H e l p ; 10 140 500 100 General settings...; showtxt 10 10 475 300 @@ GENERAL_SETTINGS ; closedlg frame 190 10 355 120 Observed R ; showtxt 205 35 85 20 Sample up to ; spin 10 10000 sample 290 35 180 20 reconsts. for each character ; subdlg 460 55 80 20 Files... ; 60 90 340 220 Save/read reconstructions... ; choose rdsvtype 20 20 260 20 Do not use files 20 40 260 20 Save reconstructions to a file for later use 20 60 280 20 Use reconstructions previously saved to a file ; showtxt 20 85 280 20 NOTE! showtxt 20 105 300 20 Please use file names without spaces in their path!!; closedlg ; check randwts 205 70 120 20 randomize radius ; + spin 1 100 reccycles 325 70 40 20 times ; = subdlg 420 100 110 20 H e l p ; 190 10 500 100 Observed r; showtxt 10 10 475 300 @@ OBSERVED_R ; closedlg frame 190 140 120 145 Randomized R ; gettxt stgetnumrepls 220 160 50 20 ; showtxt 220 180 70 30 Number of replications; check twotailed 220 225 70 20 two-tailed ; subdlg 205 260 95 20 H e l p ; 190 140 500 100 Randomization ...; showtxt 10 10 475 300 @@ RANDOMIZATION ; closedlg frame 330 140 215 145 Type of test ; choose testdown 340 155 170 18 Does X make Y change? + check stasis 380 175 150 18 Consider stasis ; = 340 200 180 18 Does Y change only with X? ; check [('havtreetags')] userlengths 340 230 180 18 Use supplied branch lengths ; subdlg 425 260 110 20 H e l p ; 300 140 300 120 Type of test; showtxt 10 10 275 300 @@ TYPE_OF_TEST ; closedlg subdlg 40 295 190 20 Significant digits and ranges .... ; 60 90 280 140 Change to no more than... ; spin 0 3 numsigdigs 20 20 180 20 significant digits; check averagit 20 45 190 20 Average min/max for terminals; closedlg closedlg if ( 'rdsvtype' ) if ( 'rdsvtype' == 1 ) getfname savefilename write Save reconstructions to ... ; else getfname savefilename read Read reconstructions from ... ; end end if ( ( 'numsigdigs' < 3 ) || 'averagit' ) var : q r ; macfloat 'numsigdigs' ; report- ; if ( 'averagit' ) loop 0 ntax if ( contmins['userx' #1] != 65.535 ) set q contmins['userx' #1] + contmaxs['userx' #1] / 2 ; xread ='userx' #1 'q'; end if ( contmins['usery' #1] != 65.535 ) set q contmins['usery' #1] + contmaxs['usery' #1] / 2 ; xread ='usery' #1 'q'; end stop else loop 0 ntax if ( contmins['userx' #1] != 65.535 ) set q contmins['userx' #1] ; set r contmaxs['userx' #1] ; xread ='userx' #1 'q'-'r' ; end if ( contmins['usery' #1] != 65.535 ) set q contmins['usery' #1] ; set r contmaxs['usery' #1] ; xread ='usery' #1 'q'-'r' ; end stop end macfloat 6 ; end if ( exstatus ) set minincx ( $stgetminx ) ; set minincy ( $stgetminy ) ; set radius ( $stgetradius ) ; set delfac ( $stgetdelfac ) ; set numrepls ( $stgetnumrepls ) ; if ( 'minincy' == 0 ) set minincy 0.0001 ; end if ( 'minincx' == 0 ) set minincx 0.0001 ; end return 1 ; else return 0 ; end end if ( ( 'numsigdigs' < 3 ) || 'averagit' ) var : q r ; macfloat 'numsigdigs' ; report- ; if ( 'averagit' ) loop 0 ntax if ( contmins['userx' #1] != 65.535 ) set q contmins['userx' #1] + contmaxs['userx' #1] / 2 ; xread ='userx' #1 'q'; end if ( contmins['usery' #1] != 65.535 ) set q contmins['usery' #1] + contmaxs['usery' #1] / 2 ; xread ='usery' #1 'q'; end stop else loop 0 ntax if ( contmins['userx' #1] != 65.535 ) set q contmins['userx' #1] ; set r contmaxs['userx' #1] ; xread ='userx' #1 'q'-'r' ; end if ( contmins['usery' #1] != 65.535 ) set q contmins['usery' #1] ; set r contmaxs['usery' #1] ; xread ='usery' #1 'q'-'r' ; end stop end macfloat 6 ; end return 1 ; /********************************************************************\ /********************************************************************\ | | | C O R R E L A T E U P | | | | This uses the command "iterrecs" to generate reconstructions, | | and travtree to move along tree | | | \********************************************************************/ label CORRELATE_UP ; set reconum 0 ; silent = all ; var : numrecsx numrecsy numrecs daprob ; if ( 'dograph' ) var: comesfrom[ (nnodes['usertree']+1) 2 ] ; end if ( 'rdsvtype' == 2 ) proc $savefilename ; else /********* Randomly select some reconstrucions for X (for subsequent testing) ***************/ set numstordx (-1) ; iterrecs 'usertree' 'userx' +/'maxstorrays' incx set numstordx ++ ; loop 0 nnodes['usertree'] set stordx['numstordx' #1] 'incx[#1]' ; stop endrecs /*************** Likewise, for Y ***************/ set numstordy (-1) ; iterrecs 'usertree' 'usery' +/'maxstorrays' incy set numstordy ++ ; loop 0 nnodes['usertree'] set stordy['numstordy' #1] 'incy[#1]' ; stop endrecs end /******************************************************* Save to a file for later use, if so requested... */ if ( 'rdsvtype' == 1 ) var : i ; macfloat 3 ; log $savefilename ; sil - file ; quote sil - all ., if ( &39maxstorrays&39 != 'maxstorrays ) errmsg Sorry, must be using the same number of reconstructions! (you saved 'maxstorrays', not &39maxstorrays&39) &38 10., end ., ; set i ntax + 1 ; quote if ( ( ntax + 1 ) != 'i' ) errmsg Sorry, must be using the same number of taxa! (you saved reconstructions for 'i') &38 10., end ., ; quote &13set numstordx 'numstordx' ., ; quote &13set numstordy 'numstordy' ., ; var stordx * ; var stordy * ; quote sil = all., proc/., ; log / ; macfloat 6 ; end /************************************************************** Note: if not using a random radius, a constant value of 'tmpradius' will be used in all cases */ set tmpradius 'radius' ; /************************************************************** Now, begin enumerating combinations of reconsts. */ set numrecsx ('numstordx'+1) ; set numrecsy ('numstordy'+1) ; set numrecs 'numrecsx' * 'numrecsy' ; sil-all ; quote TOTAL RECONSTRUCTIONS TO LOOK AT: 'numrecs' ('numrecsx'x'numrecsy')&10; sil = all ; set lor 100 ; set hir (-100) ; loop =recx 0 'numstordx' loop =noxnum 0 nnodes['usertree'] set incx[ #noxnum ] 'stordx[ #recx #noxnum ]' ; stop loop =recy 0 'numstordy' loop =noynum 0 nnodes['usertree'] set incy[ #noynum ] 'stordy[ #recy #noynum ]' ; stop set reconum ++ ; if ( 'dograph' ) if ( 'reconum' != 'dograph' ) proc/; end end progress 'reconum' 'numrecs' Checking... ; loop =cycles 1 'reccycles' set numpairs (-1) ; loop =node 0 nnodes['usertree'] set visitd[#node] 0 ; stop travtree up terms 'usertree' root - outnod if ( 'visitd['outnod']' ) proc/; end set myx 'incx['outnod']'; if ( 'randwts' ) set tmpradius getrandom[ 0 'radius' ] ; end if ( ( 'myx' >= 'minincx' ) || ( 'myx' <= (-'minincx') ) ) set didpair 0 ; set distosrc[anc['usertree' 'outnod']] (-'blength['outnod']') ; set drag[anc['usertree' 'outnod']] 0 ; travtree up terms 'usertree' 'outnod' innod set distosrc['innod'] 'distosrc[anc['usertree' 'innod']]' + 'blength['innod']' ; set day 'incy['innod']' ; set drag['innod'] 'incx['innod']' + 'drag[anc['usertree' 'innod']]' ; if ( ( 'day' >= 'minincy' ) || ( 'day' <= (-'minincy' ) ) || ( 'distosrc['innod']' > 'tmpradius' ) || ( 'innod' < ntax ) ) set didpair 1 ; set numpairs ++ ; set xpairs['numpairs'] 'drag['innod']' ; set ypairs['numpairs'] 'day' ; if ( 'dograph' ) set comesfrom['numpairs' 0] 'outnod' ; set comesfrom['numpairs' 1] 'innod' ; end set wts['numpairs'] 1 / ( 1 + ( 'distosrc['innod']' * 'delfac' ) ) ; travtree path 'usertree' 'innod' 'outnod' i set visitd['i'] 1 ; endtrav skipdes ; end endtrav if ( !'didpair' ) set numpairs ++ ; set xpairs['numpairs'] 'myx' ; set ypairs['numpairs'] 'incy['outnod']' ; set wts['numpairs'] 1 ; if ( 'dograph' ) set comesfrom['numpairs' 0] 'outnod' ; set comesfrom['numpairs' 1] 'innod' ; end proc/; end else if ( !'stasis' ) proc/; end set distosrc[anc['usertree' 'outnod']] (-'blength['outnod']') ; travtree up terms 'usertree' 'outnod' innod set dax 'incx['innod']' ; if ( ( 'dax' >= 'minincx' ) || ( 'dax' <= (-'minincx' ) ) ) skipdes ; proc/; end set distosrc['innod'] 'distosrc[anc['usertree' 'innod']]' + 'blength['innod']' ; set day 'incy['innod']' ; if ( ( 'day' >= 'minincy' ) || ( 'day' <= (-'minincy' ) ) || ( 'distosrc['innod']' > 'tmpradius' ) ) set numpairs ++ ; set xpairs['numpairs'] 'myx' ; set ypairs['numpairs'] 'day' ; if ( 'dograph' ) set comesfrom['numpairs' 0] 'outnod' ; set comesfrom['numpairs' 1] 'innod' ; end set wts['numpairs'] 1 ; travtree path 'usertree' 'innod' 'outnod' i set visitd['i'] 1 ; endtrav skipdes ; end endtrav end endtrav set numpairs ++ ; if ( 'numpairs' > 3 ) if ( 'dograph' ) /****** Graphing ... *****/ sil - all ; quote &10HAVE 'numpairs' PAIRS&10; nak-; tp; nak=; ttag=; tp; var: daan dades ; loop =node 0 nnodes['usertree'] set visitd[#node] 0 ; stop loop =pair 0 ('numpairs'-1) set daan 'comesfrom[#pair 0]' ; set dades 'comesfrom[#pair 1]' ; if ( 'dades' != 'daan' ) ttag +'dades' <#pair:'/+.3incx['dades']','/+.3incy['dades']'('/.3wts[#pair]'); travtree path 'usertree' 'dades' 'daan' - i if ( ( !'visitd['i']' ) && ( 'incx['i']' != 0 ) ) ttag +'i' '/+.3incx['i']'; set visitd['i'] 1 ; end if ( 'i' == 'daan' ) ttag +'i' |; else ttag +'i' <; end endtrav else ttag +'dades' X#pair:'/+.3incx[#pair]'/'/+.3incy[#pair]'('/.3wts[#pair]'); end stop ttag; ttag - ; loop =pair 0 ('numpairs'-1) quote Pair #pair ('/3.0comesfrom[#pair 0]','/3.0comesfrom[#pair 1]'): '/+.3xpairs[#pair]','/+.3ypairs[#pair]'&10; stop quote &10; var: tmp['numpairs']; loop =pair 0 ('numpairs'-1) set tmp[#pair] 'xpairs[#pair]' ; stop var + 60 25 tmp ypairs ; end var & xpairs ypairs *wts 'numpairs' ; if ( 'dograph' ) killrecs ; end if ( regr < 'lor' ) set lor regr ; set reclo 'reconum' ; end if ( regr > 'hir' ) set hir regr ; set rechi 'reconum' ; end end stop stop stop progress/; if ( !'dograph' ) if ( 'lor' > 10 ) /* in this case, never did any pairs!! */ sil - all ; set lor 0 ; set hir 0 ; quote CANNOT ESTABLISH CORRELATION.&10(TOO FEW PAIRS)&10; proc/; end if ( ( 'lor' <= 0 ) && ( 'hir' >= 0 ) ) sil - all ; quote RANGE OF OBSERVED R VALUES INCLUDES ZERO! ; proc/; end if ( 'lor' < 0 ) set worsr 'hir' ; else set worsr 'lor' ; end silent - all ; quote OBSERVED RANGE OF R: '/.4lor' - '/.4hir' (recs. 'reclo' - 'rechi')&10 ; sil = all ; end proc/; /**********************************************************\ /**********************************************************\ | | | R A N D O M I Z E U P | | | \**********************************************************/ label RANDOMIZE_UP ; set reconum 0 ; silent = all ; var : lapstime nocalcs i j numrecs daprob xpick ypick rali[ (nnodes['usertree']+1)] ; set nocalcs 0 ; sil-all ; quote RANDOMIZING RECONSTRUCTIONS: &10; sil = all ; set matches 1 ; set logged 1 ; set tmpradius 'radius' ; loop 1 'numrepls' set numpairs (-1) ; loop 0 nnodes['usertree'] set visitd[#2] 0 ; stop /***************************** Pick a reconstruction of X, at random, out of the 'numstordx' reconstructions stored, and copy it to the 'incx'. Also, pick a reconstruction of Y, at random, out of the 'numstordy' reconstructions stored, and randomly redistribute it over the tree (using rali, a random list of the branches) ***/ set rali randomlist[ (nnodes['usertree']+1) ] ; set xpick getrandom[ 0 'numstordx'] ; set ypick getrandom[ 0 'numstordy'] ; loop 0 nnodes['usertree'] set incx[#2] 'stordx['xpick' #2]' ; set incy['rali[#2]'] 'stordy['ypick' #2]' ; stop /************************************************************ Now, establish correlation for those reconstructions, just like done for the observed correlation, but using the randomized Y-increments */ travtree up terms 'usertree' root - outnod if ( 'visitd['outnod']' ) proc/; end set myx 'incx['outnod']'; if ( 'randwts' ) set tmpradius getrandom[ 0 'radius'] ; end if ( ( 'myx' >= 'minincx' ) || ( 'myx' <= (-'minincx') ) ) set didpair 0 ; set distosrc[anc['usertree' 'outnod']] (-'blength['outnod']') ; set drag[anc['usertree' 'outnod']] 0 ; travtree up terms 'usertree' 'outnod' innod set distosrc['innod'] 'distosrc[anc['usertree' 'innod']]' + 'blength['innod']' ; set day 'incy['innod']' ; set drag['innod'] 'incx['innod']' + 'drag[anc['usertree' 'innod']]' ; if ( ( 'day' >= 'minincy' ) || ( 'day' <= (-'minincy' ) ) || ( 'distosrc['innod']' > 'tmpradius' ) || ( 'innod' < ntax ) ) set didpair 1 ; set numpairs ++ ; set xpairs['numpairs'] 'drag['innod']' ; set ypairs['numpairs'] 'day' ; set wts['numpairs'] 1 / ( 1 + ( 'distosrc['innod']' * 'delfac' ) ) ; travtree path 'usertree' 'innod' 'outnod' i set visitd['i'] 1 ; endtrav skipdes ; end endtrav if ( !'didpair' ) set numpairs ++ ; set xpairs['numpairs'] 'myx' ; set ypairs['numpairs'] 'incy['outnod']' ; set wts['numpairs'] 1 ; proc/; end else if ( !'stasis' ) proc/; end set distosrc[anc['usertree' 'outnod']] (-'blength['outnod']') ; travtree up terms 'usertree' 'outnod' innod set dax 'incx['innod']' ; if ( ( 'dax' >= 'minincx' ) || ( 'dax' <= (-'minincx' ) ) ) skipdes ; proc/; end set distosrc['innod'] 'distosrc[anc['usertree' 'innod']]' + 'blength['innod']' ; set day 'incy['innod']' ; if ( ( 'day' >= 'minincy' ) || ( 'day' <= (-'minincy' ) ) || ( 'distosrc['innod']' > 'tmpradius' ) ) set numpairs ++ ; set xpairs['numpairs'] 'myx' ; set ypairs['numpairs'] 'day' ; set wts['numpairs'] 1 ; travtree path 'usertree' 'innod' 'outnod' i set visitd['i'] 1 ; endtrav skipdes ; end endtrav end endtrav set numpairs ++ ; if ( 'numpairs' > 3 ) var & xpairs ypairs *wts 'numpairs' ; /******************************************************* Establish regression for this randomization, and compare with observed value of r */ if ( 'twotailed' ) if ( regr < 0 ) if ( ( 'ispos' && ( (-regr) >= 'lor' ) ) || ( (!'ispos') && ( regr <= 'hir' ) ) ) set matches ++ ; end else if ( ( (!'ispos') && ( (-regr) <= 'hir' ) ) || ( 'ispos' && ( regr >= 'lor' ) ) ) set matches ++ ; end end else if ( 'ispos' && ( regr >= 'lor' ) ) set matches ++ ; else if ( (!'ispos') && ( regr <= 'hir' ) ) set matches ++ ; end end end set signif ( 'matches' / ( #1 + 1 ) ) ; set i #1 ; set j regr ; set lapstime time ; if ( 'logged' == 10 ) sil - all ; else if ( 'lapstime' >= 1 ) resettime ; sil - con ; end end if ( !windows ) quote &0 '/5.0i' REPLS. - P(inst ) '/.3signif' ('matches' matches) - this R = '/+j'; if ( 'logged' == 10 ) quote &10 ; set logged 0 ; end else progress 'i' 'numrepls' 'i' REPLS - P(inst) '/.3signif' ('matches' matches) ; if ( 'logged' == 10 ) quote &32 &32 &32 '/5.0i' REPLS. - P(inst ) '/.3signif' ('matches' matches) - this R = '/+j'&10; set logged 0 ; end end set logged ++ ; sil = all ; else set nocalcs ++ ; setloop #1 ; end stop if ( windows ) progress/; end lquote - ; sil - all ; if ( !windows ) quote &0 '/5.0i' REPLS. - P(final) '/.3signif' ('matches' matches) &10; else quote &32 &32 &32 '/5.0i' REPLS. - P(final) '/.3signif' ('matches' matches) &10; end proc/; /****************************************************************\ /****************************************************************\ | | | C O R R E L A T E D O W N | | | \****************************************************************/ label CORRELATE_DOWN ; set reconum 0 ; silent = all ; var : distoy numrecsx numrecsy foundanx numrecs daprob ; if ( 'dograph' ) var: comesfrom[ (nnodes['usertree']+1) 2 ] ; end if ( 'rdsvtype == 2 ) proc $savefilename ; else /********* Randomly select some reconstrucions for X (for subsequent testing) ***************/ set numstordx (-1) ; iterrecs 'usertree' 'userx' +/'maxstorrays' incx set numstordx ++ ; loop 0 nnodes['usertree'] set stordx['numstordx' #1] 'incx[#1]' ; stop endrecs /*************** Likewise, for Y ***************/ set numstordy (-1) ; iterrecs 'usertree' 'usery' +/'maxstorrays' incy set numstordy ++ ; loop 0 nnodes['usertree'] set stordy['numstordy' #1] 'incy[#1]' ; stop endrecs end /******************************************************* Save to a file for later use, if so requested... */ if ( 'rdsvtype' == 1 ) var : i ; macfloat 3 ; log $savefilename ; sil - file ; quote sil - all ., if ( &39maxstorrays&39 != 'maxstorrays ) errmsg Sorry, must be using the same number of reconstructions! (you saved 'maxstorrays', not &39maxstorrays&39) &38 10., end ., ; set i ntax + 1 ; quote if ( ( ntax + 1 ) != 'i' ) errmsg Sorry, must be using the same number of taxa! (you saved reconstructions for 'i') &38 10., end ., ; quote &13set numstordx 'numstordx' ., ; quote &13set numstordy 'numstordy' ., ; var stordx * ; var stordy * ; quote sil = all., proc/., ; log / ; macfloat 6 ; end /************************************************************** Note: if not using a random radius, a constant value of 'tmpradius' will be used in all cases */ set tmpradius 'radius' ; /************************************************************** Now, begin enumerating combinations of reconsts. */ set numrecsx ('numstordx'+1) ; set numrecsy ('numstordy'+1) ; set numrecs 'numrecsx' * 'numrecsy' ; sil-all ; quote TOTAL RECONSTRUCTIONS TO LOOK AT: 'numrecs' ('numrecsx'x'numrecsy')&10; sil = all ; set lor 100 ; set hir (-100) ; loop =recx 0 'numstordx' loop =nonum 0 nnodes['usertree'] set incx[ #nonum ] 'stordx[ #recx #nonum ]' ; stop loop =recy 0 'numstordy' loop =nonum 0 nnodes['usertree'] set incy[ #nonum ] 'stordy[ #recy #nonum ]' ; stop set reconum ++ ; if ( 'dograph' ) if ( 'reconum' != 'dograph' ) proc/; end end progress 'reconum' 'numrecs' Checking... ; loop =cycle 1 'reccycles' set numpairs (-1) ; set tmpradius 'radius' ; loop =node 0 nnodes['usertree'] if ( #node == root ) continue ; end if ( 'randwts' ) set tmpradius getrandom [ 0 'radius' ] ; end set day 'incy[ #node ]' ; set outnod #node ; if ( ( 'day' >= 'minincy' ) || ( 'day' <= (-'minincy' ) ) ) set i 0 ; set dax 0 ; set foundanx 0 ; set distoy 0 ; travtree path 'usertree' 'outnod' root innod set myx 'incx['innod']'; set dax += 'myx' ; if ( !'foundanx' ) if ( ( 'myx' >= 'minincx' ) || ( 'myx' <= (-'minincx' ) ) ) set distoy 'i' ; set foundanx 1 ; end end set i ++ ; if ( 'i' > 'tmpradius' ) killtrav ; end endtrav set numpairs ++ ; set xpairs [ 'numpairs' ] 'dax' ; set ypairs [ 'numpairs' ] 'day' ; if ( 'dograph' ) set comesfrom[ 'numpairs' 0 ] 'innod' ; set comesfrom[ 'numpairs' 1 ] 'outnod' ; end if ( 'foundanx' ) set wts['numpairs'] 1 / ( 1 + ( 'distoy' * 'delfac' ) ) ; else set wts['numpairs'] 1 ; end end stop set numpairs ++ ; if ( 'numpairs' > 3 ) if ( 'dograph' ) /*** Graphing ... **/ sil - all ; quote &10HAVE 'numpairs' PAIRS&10; loop =pair 0 ('numpairs'-1) quote Pair #pair ('/3.0comesfrom[#pair 0]','/3.0comesfrom[#pair 1]'): '/+.3xpairs[#pair]','/+.3ypairs[#pair]'&10; stop quote &10; var: tmp['numpairs']; loop =pair 0 ('numpairs'-1) set tmp[#pair] 'xpairs[#pair]' ; stop var + 60 25 tmp ypairs ; end var & xpairs ypairs *wts 'numpairs' ; if ( 'dograph' ) killrecs ; end if ( regr < 'lor' ) set lor regr ; set reclo 'reconum' ; end if ( regr > 'hir' ) set hir regr ; set rechi 'reconum' ; end end stop stop stop progress/; if ( !'dograph' ) if ( 'lor' > 10 ) /* in this case, never did any pairs!! */ sil - all ; set lor 0 ; set hir 0 ; quote CANNOT ESTABLISH CORRELATION.&10(TOO FEW PAIRS)&10; proc/; end if ( ( 'lor' <= 0 ) && ( 'hir' >= 0 ) ) sil - all ; quote RANGE OF OBSERVED R VALUES INCLUDES ZERO! ; proc/; end if ( 'lor' < 0 ) set worsr 'hir' ; else set worsr 'lor' ; end silent - all ; quote OBSERVED RANGE OF R: '/.4lor' - '/.4hir' (recs. 'reclo' - 'rechi')&10 ; sil = all ; end proc/; /**************************************************************\ /**************************************************************\ | | | R A N D O M I Z E D O W N | | | \**************************************************************/ label RANDOMIZE_DOWN ; set reconum 0 ; silent = all ; var : lapstime nocalcs i j numrecs daprob xpick ypick rali[ (nnodes['usertree']+1)] foundanx distoy ; set nocalcs 0 ; sil-all ; quote RANDOMIZING RECONSTRUCTIONS: &10; sil = all ; set matches 1 ; set logged 1 ; set tmpradius 'radius' ; loop 1 'numrepls' set numpairs (-1) ; set tmpradius 'radius' ; set rali randomlist[ (nnodes['usertree']+1) ] ; set xpick getrandom[ 0 'numstordx'] ; set ypick getrandom[ 0 'numstordy'] ; loop 0 nnodes['usertree'] set incx[#2] 'stordx['xpick' #2]' ; set incy['rali[#2]'] 'stordy['ypick' #2]' ; stop loop =node 0 nnodes['usertree'] if ( #node == root ) continue ; end if ( 'randwts' ) set tmpradius getrandom [ 0 'radius' ] ; end set day 'incy[ #node ]' ; set outnod #node ; if ( ( 'day' >= 'minincy' ) || ( 'day' <= (-'minincy' ) ) ) set i 0 ; set dax 0 ; set foundanx 0 ; set distoy 0 ; travtree path 'usertree' 'outnod' root innod set myx 'incx['innod']'; set dax += 'myx' ; if ( !'foundanx' ) if ( ( 'myx' >= 'minincx' ) || ( 'myx' <= (-'minincx' ) ) ) set distoy 'i' ; set foundanx 1 ; end end set i ++ ; if ( 'i' > 'tmpradius' ) killtrav ; end endtrav set numpairs ++ ; set xpairs [ 'numpairs' ] 'dax' ; set ypairs [ 'numpairs' ] 'day' ; if ( 'dograph' ) set comesfrom[ 'numpairs' 0 ] 'innod' ; set comesfrom[ 'numpairs' 1 ] 'outnod' ; end if ( 'foundanx' ) set wts['numpairs'] 1 / ( 1 + ( 'distoy' * 'delfac' ) ) ; else set wts['numpairs'] 1 ; end end stop set numpairs ++ ; if ( 'numpairs' > 3 ) var & xpairs ypairs *wts 'numpairs' ; if ( 'twotailed' ) if ( regr < 0 ) if ( ( 'ispos' && ( (-regr) >= 'lor' ) ) || ( (!'ispos') && ( regr <= 'hir' ) ) ) set matches ++ ; end else if ( ( (!'ispos') && ( (-regr) <= 'hir' ) ) || ( 'ispos' && ( regr >= 'lor' ) ) ) set matches ++ ; end end else if ( 'ispos' && ( regr >= 'lor' ) ) set matches ++ ; else if ( (!'ispos') && ( regr <= 'hir' ) ) set matches ++ ; end end end set signif ( 'matches' / ( #1 + 1 ) ) ; set i #1 ; set j regr ; set lapstime time ; if ( 'logged' == 10 ) sil - all ; else if ( 'lapstime' >= 1 ) resettime ; sil - con ; end end if ( !windows ) quote &0 '/5.0i' REPLS. - P(inst ) '/.3signif' ('matches' matches) - this R = '/+j'; if ( 'logged' == 10 ) quote &10 ; set logged 0 ; end else progress 'i' 'numrepls' 'i' REPLS - P(inst) '/.3signif' ('matches' matches) ; if ( 'logged' == 10 ) quote &32 &32 &32 '/5.0i' REPLS. - P(inst ) '/.3signif' ('matches' matches) - this R = '/+j'&10; set logged 0 ; end end set logged ++ ; sil = all ; else set nocalcs ++ ; setloop #1 ; end stop if ( windows ) progress/; end lquote - ; sil - all ; if ( !windows ) quote &0 '/5.0i' REPLS. - P(final) '/.3signif' ('matches' matches) &10; else quote &32 &32 &32 '/5.0i' REPLS. - P(final) '/.3signif' ('matches' matches) &10; end proc/; /************************************************************\ /************************************************************\ | | | C H E C K P A R M S | | (it makes sure parameters are ok) | | | \************************************************************/ label CHECKPARMS ; if ( (!iscont['userx']) || (!iscont['usery']) ) errmsg Sorry, both characters X and Y must be continuous ; end if ( tsize['usertree'] != root ) errmsg Sorry, all taxa in the matrix must be included in the tree!!; end proc/; /***************************************************************\ /***************************************************************\ | | | D I S P L A Y H E L P | | | \***************************************************************/ label DISPLAY_HELP lquote [ ; sil - all ; quote * This script tests character correlation between continuous characters. Argument: default: tree N 0 chars X Y 0,1 minincx N 0 minincy N 0 radius N 4 stasis no delfac N 1 sample N 50 randwts no cycles N 1, if randwts: 5 repls N 100 userlengths no twotailed no testdown up dograph N no notest test average no sigdigs N 3 saveto filename readfrom filename &10; proc/; /***************************************************************\ /***************************************************************\ | | | T E X T F O R D I A L O G S | | | \***************************************************************/ label CHOOSE_CHARS Note that the tree must be complete (i.e. have all taxa). Both X and Y must be continuous characters ; @@ label GENERAL_SETTINGS Changes below Min/Max increases are ignored. Max radius is the number of nodes beyond which the dependent variable is considered to not have changed. As the delay factor increases, changes in the dependent variable occuring more nodes away are considered to be less compelling evidence of correlation. @@ label OBSERVED_R The observed R is calculated by using a sample of reconstructions for the dependent and independent variables (the lowest R is used). The radius can be randomized for each combination of reconstructions, so that the conclusion does not depend on a specific radius value. @@ label RANDOMIZATION The randomization re-distributes observed changes (for different alternative reconstructions for X and Y) on the tree, and recalculates R. @@ label TYPE_OF_TEST Two types of test can be performed. Branch-lengths must be provided as tree-tags. @@