sil=all ; screen 39x80 ; if ( ntrees < 0 ) errmsg Must read data (and trees) before running ; end macro *10 ( ( ( 2*ntax+2 ) * 35 ) + 1000 ) ; macro [ ( ( ( ( ( 2*ntax+2 ) * 35 ) + 1000 ) / 1000 ) + 1000 ) ; macro=; /*** Please note !! There is an update for this script: http://www.zmuc.dk/public/phylogeny/TNT/scripts/xdelcor.run That script is better at handling large data sets with very numerous reconstructions (and can also save the arrays of increments/decrements to a file, so that testing with different options can be done without having to generate again all possible reconstructions for each of the characters). It requires a version of TNT for March, 2010, or later. The present script is --this note aside-- exactly the same one as used by Giannini & Goloboff (2010), and is preserved as such 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 | | calls (with "goto") the rest of the routines | | that are used in the script. | | | \************************************************************/ var: 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 ; set maxstorrays 10 ; /* Change this if you want to store more arrays! */ var: stordx [ 'maxstorrays' (2*ntax+2) ] stordy [ 'maxstorrays' (2*ntax+2) ] ; /******************** 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 ; 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 100 ; /* 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 */ 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 errmsg Wrong argument: %('i'); 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 260 120 Observed R ; spin 10 10000 sample 205 30 150 20 Sample of reconstruction ; showtxt 290 47 180 20 combinations ; check randwts 205 70 120 20 randomize radius ; + spin 1 100 reccycles 325 70 40 20 times ; = subdlg 330 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 closedlg 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 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 set numrecsx 0 ; set numrecsy 0 ; iterrecs 'usertree' 'userx' incx set numrecsx ++ ; endrecs iterrecs 'usertree' 'usery' incy set numrecsy ++ ; endrecs set numrecs 'numrecsx' * 'numrecsy' ; sil-all ; quote TOTAL RECONSTRUCTIONS: 'numrecs' ('numrecsx'x'numrecsy')&10; sil = all ; /******************************************************************** Calculate the probability P of rejecting a combination, such that P is expected to produce the desired number of combinations (='sample'). Note this can't be less than 1/1000) */ if ( 'sample' < 10 ) set sample 10 ; end set daprob ( 'sample' / 'numrecs' ) * 1000 ; if ( 'daprob' < 1 ) set daprob 1 ; end if ( 'numrecs' < 20 ) set daprob 1000 ; end set lor 100 ; set hir (-100) ; /********* Randomly select some reconstrucions for X (for subsequent testing) ***************/ set numstordx (-1) ; if ( 'numrecsx' <= 'maxstorrays' ) iterrecs 'usertree' 'userx' + incx set numstordx ++ ; loop 0 nnodes['usertree'] set stordx['numstordx' #1] 'incx[#1]' ; stop endrecs else var : atmpsel['maxstorrays'] ; loop 1 'maxstorrays' set i getrandom[0 ('numrecsx'-1)] ; if ( 'numstordx' >= 0 ) set matches 0 ; loop 0 'numstordx' if ( 'atmpsel[#2]' == 'i' ) set matches 1 ; endloop ; end stop if ( 'matches' ) setloop #1 ; end end set numstordx ++ ; set atmpsel['numstordx'] 'i' ; stop set numstordx (-1) ; set i 0 ; iterrecs 'usertree' 'userx' + incx set matches 0 ; loop 0 ('maxstorrays'-1) if ( 'atmpsel[#1]' == 'i' ) set matches 1 ; endloop ; end stop set i ++ ; if ( !'matches' ) proc/; end set numstordx ++ ; loop 0 nnodes['usertree'] set stordx['numstordx' #1] 'incx[#1]' ; stop endrecs var - atmpsel ; end /*************** Likweise, for Y ***************/ set numstordy (-1) ; if ( 'numrecsy' <= 'maxstorrays' ) iterrecs 'usertree' 'usery' + incy set numstordy ++ ; loop 0 nnodes['usertree'] set stordy['numstordy' #1] 'incy[#1]' ; stop endrecs else var : atmpsel['maxstorrays'] ; loop 1 'maxstorrays' set i getrandom[0 ('numrecsy'-1)] ; if ( 'numstordy' >= 0 ) set matches 0 ; loop 0 'numstordy' if ( 'atmpsel[#2]' == 'i' ) set matches 1 ; endloop ; end stop if ( 'matches' ) setloop #1 ; end end set numstordy ++ ; set atmpsel['numstordy'] 'i' ; stop set numstordy (-1) ; set i 0 ; iterrecs 'usertree' 'usery' + incy set matches 0 ; loop 0 ('maxstorrays'-1) if ( 'atmpsel[#1]' == 'i' ) set matches 1 ; endloop ; end stop set i ++ ; if ( !'matches' ) proc/; end set numstordy ++ ; loop 0 nnodes['usertree'] set stordy['numstordy' #1] 'incy[#1]' ; stop endrecs var - atmpsel ; 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. */ iterrecs 'usertree' 'userx' + incx iterrecs 'usertree' 'usery' + incy set reconum ++ ; if ( !'dograph' ) if ( getrandom[ 1 1000 ] > 'daprob' ) proc/; end else 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 endrecs endrecs 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 set numrecsx 0 ; set numrecsy 0 ; iterrecs 'usertree' 'userx' incx set numrecsx ++ ; endrecs iterrecs 'usertree' 'usery' incy set numrecsy ++ ; endrecs set numrecs 'numrecsx' * 'numrecsy' ; sil-all ; quote TOTAL RECONSTRUCTIONS: 'numrecs' ('numrecsx'x'numrecsy')&10; sil = all ; if ( 'sample' < 10 ) set sample 10 ; end set daprob ( 'sample' / 'numrecs' ) * 1000 ; if ( 'daprob' < 1 ) set daprob 1 ; end if ( 'numrecs' < 20 ) set daprob 1000 ; end set lor 100 ; set hir (-100) ; /***** Randomly select some reconstrucions for X (for subsequent testing) *******/ set numstordx (-1) ; if ( 'numrecsx' <= 'maxstorrays' ) iterrecs 'usertree' 'userx' + incx set numstordx ++ ; loop 0 nnodes['usertree'] set stordx['numstordx' #1] 'incx[#1]' ; stop endrecs else var : atmpsel['maxstorrays'] ; loop 1 'maxstorrays' set i getrandom[0 ('numrecsx'-1)] ; if ( 'numstordx' >= 0 ) set matches 0 ; loop 0 'numstordx' if ( 'atmpsel[#2]' == 'i' ) set matches 1 ; endloop ; end stop if ( 'matches' ) setloop #1 ; end end set numstordx ++ ; set atmpsel['numstordx'] 'i' ; stop set numstordx (-1) ; set i 0 ; iterrecs 'usertree' 'userx' + incx set matches 0 ; loop 0 ('maxstorrays'-1) if ( 'atmpsel[#1]' == 'i' ) set matches 1 ; endloop ; end stop set i ++ ; if ( !'matches' ) proc/; end set numstordx ++ ; loop 0 nnodes['usertree'] set stordx['numstordx' #1] 'incx[#1]' ; stop endrecs var - atmpsel ; end /***** Likweise, for Y *************/ set numstordy (-1) ; if ( 'numrecsy' <= 'maxstorrays' ) iterrecs 'usertree' 'usery' + incy set numstordy ++ ; loop 0 nnodes['usertree'] set stordy['numstordy' #1] 'incy[#1]' ; stop endrecs else var : atmpsel['maxstorrays'] ; loop 1 'maxstorrays' set i getrandom[0 ('numrecsy'-1)] ; if ( 'numstordy' >= 0 ) set matches 0 ; loop 0 'numstordy' if ( 'atmpsel[#2]' == 'i' ) set matches 1 ; endloop ; end stop if ( 'matches' ) setloop #1 ; end end set numstordy ++ ; set atmpsel['numstordy'] 'i' ; stop set numstordy (-1) ; set i 0 ; iterrecs 'usertree' 'usery' + incy set matches 0 ; loop 0 ('maxstorrays'-1) if ( 'atmpsel[#1]' == 'i' ) set matches 1 ; endloop ; end stop set i ++ ; if ( !'matches' ) proc/; end set numstordy ++ ; loop 0 nnodes['usertree'] set stordy['numstordy' #1] 'incy[#1]' ; stop endrecs var - atmpsel ; end /*** Begin enumerating combinations of reconsts. ******/ iterrecs 'usertree' 'userx' + incx iterrecs 'usertree' 'usery' + incy set reconum ++ ; if ( !'dograph' ) if ( getrandom[ 1 1000 ] > 'daprob' ) proc/; end else 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 endrecs endrecs 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 100 randwts no cycles N 1, if randwts: 5 repls N 100 userlengths no twotailed no testdown up dograph N no notest &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. @@