How to use the DDsim package

Copyright  A.Grafen and M. Ridley   1997

Introduction

This file is intended to help you use DDsim, a Mathematica package that generates simulated discrete comparative data.  For an explanation of the simulations, and the reasons behind them, see the original paper by Grafen and Ridley (1997. A new model for discrete character evolution. 
J. theor. Biol. 184, 7-14).

There are two main parts.  First, a demonstration of how to use the package.  Second, an explanation of how phylogenies have to be encoded for use with the program.

You can execute the input cells in this notebook, and use the package as described.  So this is an active notebook.  You can  evaluate each input cell as you come to it (recognising them by the lack of an extra bar on the innermost bracket over to the right of the window).  Or you can select 'Evaluate Notebook' from the menus.

This package was written to create datasets against which to test different discrete comparative methods.  The programming is rough and ready, and would win no prizes for elegance (or, probably, efficiency).  However, the algorithms do work and have been extensively tested.

Part 1 - A worked example of the use of DDsim

Loading the package

The first task is to load the package.  You can use the Needs command, as follows.  If Mathematica  can't find the file, it should prompt you to find it.  You will then have to remember which folder or directory you put it in.

Needs["DDsim`"]

Preparing to use the package

Specifying a simulation requires specifying a lot of detail.  The details are given to the command "trial", which then produces the data.  You need to call trial with the following arguments

 1. name of output file
 2. a list with ten elements
    (i) number of simulated datasets required
   (ii) initial states for C and D, as {c0,d0}
  (iii) the value of rho used to transform node heights
   (iv) the rate of evolution
    (v) the phylogeny
   (vi) rule for determining node heights
  (vii) transition-rate matrix for C
 (viii) transition-rate matrix for D
   (ix) rule for determining A from B and C
    (x) rule for determining B from D

It is likely you will want to ask for similar simulations a number of times.  If you want to know how long it will take to create 1000 datasets, it is worthwhile creating 2 and 12, to work out how long it takes to make 10 once the startup time has finished.  If you want to get a 'reasonable' amount of variation and phylogenetic association, you might want to vary rho and rate experimentally. 

So we will begin by specifying the phylogeny, the node-height rule, the transition matrices and the A and B rules.  Then we will build them into a convenient way to call for datasets with those elements fixed, while allowing the others to be specified each time.

2(v) The tree

First, the tree.  To specify your own tree, see the second part of this file.  For the moment, we will use a command that automatically creates a phylogeny in the appropriate format.  'symtree' makes a symmetrical tree.  Let's make a small tree with only 16 species, to save time.  Then the command

mytree=symtree[4,2]

makes a tree with four branches per higher node, and two levels below the root.  Incidentally, you can remind yourself of the function of commands provided in the usual Mathematica way:

?symtree

2(vi) The node height rule

Our next piece of information is the 'node-height rule'.  The branch lengths are calculated by assigning a height to each node, and then the length of any branch segment is the difference between the heights of its upper and lower nodes.  We specify a height function that takes the number of species below the node as argument, and delivers a real number as result.  If you don't understand the functional notation, either just copy blindly from this notebook, or read the Mathematica documentation about 'Function'.  A log rule would be produced by

  myhtf=N[Log[#]]&
 
Instead, we will use the 'Figure 2' method of Grafen (1989 - see Grafen and Ridley (loc. cit.) for the full reference):

myhtf=N[#-1]&

The reason for the use of 'N' is that Mathematica is capable of trying to do calculations with the numbers as rationals instead of reals.  The computations cannot normally be completed with rationals, and Mathematica will eventually seize up, going over memory or exceeding your patience.  Using 'N' tells Mathematica to take the result as a real number and not a rational.

2(vii) and 2(viii)  The transition-rate matrices

Now we specify the transition matrices for characters C and D.  See Grafen and Ridley for an explanation of what the numbers mean.  For the present, it suffices to say that the matrices must be square, and the transition matrix for C must have the same number of rows as there are states C, and similarly for D.  In our case, we will have 6 states of C and 6 states of D, and we specify

myctrans = {{-0.01, 0.01, 0., 0., 0., 0.}, {0.2, -0.82, 0.6, 0.02, 0., 0.}, 
   {0., 0.6, -0.6199999999999999999, 0.02, 0., 0.}, 
   {0., 0., 0.02, -0.6199999999999999999, 0.6, 0.}, 
   {0., 0., 0.02, 0.6, -0.82, 0.2}, {0., 0., 0., 0., 0.01, -0.01}};
 
mydtrans = {{-0.01, 0.01, 0., 0., 0., 0.}, {0.2, -0.82, 0.6, 0.02, 0., 0.}, 
   {0., 0.6, -0.6199999999999999999, 0.02, 0., 0.}, 
   {0., 0., 0.02, -0.6199999999999999999, 0.6, 0.}, 
   {0., 0., 0.02, 0.6, -0.82, 0.2}, {0., 0., 0., 0., 0.01, -0.01}};

The rates must sum to zero, so if we add up each column, it should come to zero.  We can check thus (and should get {0,0,0,0,0,0} or very close in each case):

Plus@@Transpose[myctrans]

Plus@@Transpose[mydtrans]

2(ix) and 2(x)  The A and B rules

The arule and brule are last.  The brule specifies B from D.  The form given here is simply a list.  The list has the same number of elements as there are states of D.  Each element specifies what value of B you wish to give to each element of D, in turn.  Again, the functional notation is used, but you can edit just the list and leave the rest alone.

mybrule=({1,1,2,1,2,2}[[#]])&

This rule says D=1, 2 or 4 makes B=1; while D=3, 5 or 6 makes B=2.  Next, the arule must specify A as a function of B and C.  The second argument is C.  To make data under the null hypothesis that A and B are unrelated, we simply omit any mention of the first argument (which is B) from the rule for making A, thus:

myarule = ({1, 1, 2, 1, 2, 2}[[#2]])&

Had we wanted to make non-null data, we could do this by giving a two-dimensional list like:

 myarule =
 ({{1, 1, 1, 1, 1, 2},{1, 1, 2, 1, 2, 2}}[[#1,#2]])&

This rule would specify that when B=1, A is 1 except when C=6; but when B=2, A is 1 when C=3, 5 or 6.  Thus B affects A in the cases C=3 and C=5.  (You can copy and paste this rule and then edit the list to make your own non-null data.)

Defining a convenient shortcut

The next step is to make a command that will allow us to store up the values of all the parameters specified so far, and allow us to change only the others.  We do this as

myparams[nsims_,{catroot_,datroot_},rho_,rate_]:=
 {nsims,{catroot,datroot},rho,rate,mytree,myhtf,
  myctrans,mydtrans,myarule,mybrule}

Creating the datasets

This allows us to call the main command 'trial' more simply.  Let's ask for 20 datasets, starting with C=3, D=2 at the root, with rho=0.8 and rate=3.4.  To put the results into a file called "mysims.dat", we would do:

trial["mysims.dat",myparams[20,{3,2},0.8,3.4]]

When I did this on my own machine, I got the following result:

trial["mysims.dat",myparams[20,{3,2},0.8,3.4]]

This shows that Mathematica took about 50 seconds, there is one "Null" for each dataset, and confirms that the output went into the file "mysims.dat"

Format of the output file

I then read in mysims.dat, the file created in the previous section and got this:

!!mysims.dat

The results are interpreted as follows.  The first few lines show the parameters used to create the datasets.  Then each of the twenty datasets is given, with the data for each species as {A,B}, in the same order as the species appear in the phylogeny.  The data can be read in as Mathematica expressions.  If you want to read them into some other program, then it will be necessary to edit the file in some way.

When performing simulations ourselve, Mark Ridley and I kept these output files on file and used edited versions for the simulations themselves.  The originals are important in case you need to find out what parameter values were used in their creation.

Part 2 - The coding of phylogenies

The form of phylogenies

The phylogeny supplied to 'trial' must be in a particular form.  Let's look again at mytree, as generated earlier by mytree=symtree[4,2]:

mytree

Each higher node is of the form

  {Null, lower_node_1,lower_node_2,lower_node_3,lower_node_4}
  
though of course there can be as few as 2 lower nodes and as many as you like.  The initial 'Null' is ugly but required.  Each lower node is either itself a higher node, of the same form, or a species node.  A species node is represented just as 'tip[]'.  So 'mytree' has one top node which has four lower nodes.  Each lower node is itself a higher node that contains four species.

Any rooted tree can be represented in this way.  In particular, there is no need for the depth to be the same.  Thus a ladder phylogeny with four species would look like this:

  {Null, tip[], {Null, tip[], {Null, tip[],tip[]}}}
  
If the phylogeny is wrongly formatted, then the results of applying the package are unpredictable.

The binary and insect phylogenies used by Grafen and Ridley

The binary and insect phylogenies used by Grafen and Ridley are given here, as bintree and insecttree, respectively, for reference.  Each tree has 360 species.  The other tree used in that paper was simply generated by symtree[4,4].

bintree= {Null, {Null, {Null, {Null, 
      {Null, {Null, {Null, {Null, tip[], 
          {Null, {Null, tip[], tip[]}, tip[]}}, 
         {Null, {Null, {Null, tip[], {Null, tip[], tip[]}}, tip[]}, 
          {Null, {Null, tip[], {Null, tip[], tip[]}}, tip[]}}}, 
        {Null, {Null, tip[], {Null, tip[], tip[]}}, tip[]}}, 
       {Null, {Null, {Null, {Null, {Null, tip[], tip[]}, tip[]}, 
          tip[]}, {Null, tip[], {Null, tip[], {Null, tip[], tip[]}}}}, 
        {Null, {Null, tip[], {Null, tip[], {Null, tip[], tip[]}}}, 
         {Null, tip[], {Null, tip[], {Null, tip[], tip[]}}}}}}, 
      {Null, {Null, {Null, {Null, tip[], 
          {Null, {Null, tip[], tip[]}, tip[]}}, 
         {Null, tip[], {Null, {Null, tip[], tip[]}, tip[]}}}, 
        {Null, {Null, {Null, {Null, tip[], tip[]}, tip[]}, tip[]}, 
         {Null, {Null, {Null, tip[], tip[]}, tip[]}, tip[]}}}, 
       {Null, {Null, {Null, tip[], {Null, tip[], {Null, tip[], tip[]}}}, 
         {Null, tip[], {Null, {Null, tip[], tip[]}, tip[]}}}, 
        {Null, {Null, {Null, tip[], {Null, tip[], tip[]}}, tip[]}, 
         {Null, {Null, tip[], {Null, tip[], tip[]}}, tip[]}}}}}, 
     {Null, {Null, {Null, {Null, 
         {Null, {Null, {Null, {Null, tip[], tip[]}, tip[]}, tip[]}, 
          {Null, tip[], {Null, tip[], {Null, tip[], tip[]}}}}, 
         {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}}, 
        {Null, {Null, tip[], {Null, tip[], tip[]}}, tip[]}}, 
       {Null, {Null, {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}, 
         {Null, {Null, {Null, {Null, tip[], tip[]}, 
            {Null, tip[], tip[]}}, 
           {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}}, 
          {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}}}, 
        {Null, {Null, {Null, tip[], 
           {Null, tip[], {Null, tip[], tip[]}}}, 
          {Null, {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}, 
           {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}}}, 
         {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}}}}, 
      {Null, {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}, 
       {Null, {Null, {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}, 
         {Null, tip[], {Null, {Null, tip[], tip[]}, tip[]}}}, 
        {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}}}}}, 
    {Null, {Null, {Null, {Null, 
        {Null, {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}, 
         {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}}, 
        {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}}, 
       {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}}, 
      {Null, {Null, {Null, {Null, 
          {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}, 
          {Null, {Null, tip[], {Null, tip[], tip[]}}, tip[]}}, 
         {Null, tip[], {Null, {Null, tip[], tip[]}, tip[]}}}, 
        {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}}, 
       {Null, {Null, tip[], {Null, {Null, tip[], tip[]}, tip[]}}, 
        {Null, {Null, {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}, 
          {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}}, 
         {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}}}}}, 
     {Null, {Null, {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}, 
       {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}}, 
      {Null, {Null, tip[], {Null, tip[], {Null, tip[], tip[]}}}, 
       {Null, tip[], {Null, tip[], {Null, tip[], tip[]}}}}}}}, 
   {Null, {Null, {Null, {Null, tip[], {Null, tip[], tip[]}}, tip[]}, 
     {Null, {Null, {Null, tip[], {Null, {Null, tip[], tip[]}, tip[]}}, 
       {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}}, 
      {Null, {Null, tip[], {Null, tip[], tip[]}}, tip[]}}}, 
    {Null, {Null, {Null, {Null, tip[], 
        {Null, tip[], {Null, tip[], tip[]}}}, 
       {Null, {Null, {Null, {Null, tip[], tip[]}, tip[]}, tip[]}, 
        {Null, {Null, {Null, tip[], {Null, tip[], tip[]}}, tip[]}, 
         {Null, {Null, tip[], {Null, tip[], tip[]}}, tip[]}}}}, 
      {Null, {Null, {Null, {Null, tip[], tip[]}, tip[]}, tip[]}, 
       {Null, {Null, {Null, tip[], {Null, tip[], tip[]}}, tip[]}, 
        {Null, {Null, {Null, {Null, tip[], tip[]}, tip[]}, tip[]}, 
         {Null, {Null, tip[], {Null, tip[], tip[]}}, tip[]}}}}}, 
     {Null, {Null, {Null, tip[], {Null, tip[], {Null, tip[], tip[]}}}, 
       {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}}, 
      {Null, {Null, tip[], {Null, {Null, tip[], tip[]}, tip[]}}, 
       {Null, tip[], {Null, {Null, tip[], tip[]}, tip[]}}}}}}};

insecttree={Null, {Null, tip[], tip[]}, {Null, tip[], 
   {Null, {Null, tip[], tip[], {Null, tip[], tip[]}}, 
    {Null, tip[], tip[], tip[]}}, 
   {Null, tip[], {Null, {Null, 
      {Null, {Null, {Null, {Null, 
          {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}, 
          {Null, {Null, {Null, {Null, tip[], tip[]}, tip[]}, tip[]}, 
           tip[]}}, {Null, {Null, {Null, {Null, tip[], tip[]}, tip[]}, 
           {Null, {Null, tip[], tip[]}, 
            {Null, tip[], {Null, tip[], tip[]}}}}, tip[]}}, tip[]}, tip[]}, 
      tip[]}, {Null, {Null, tip[], tip[], tip[]}, 
      {Null, {Null, {Null, {Null, tip[], tip[], 
          {Null, tip[], {Null, tip[], tip[]}}, 
          {Null, tip[], {Null, tip[], 
            {Null, {Null, {Null, {Null, tip[], tip[]}, tip[]}, tip[]}, 
             {Null, {Null, {Null, tip[], tip[], tip[]}, tip[]}, 
              {Null, {Null, {Null, {Null, tip[], tip[]}, tip[]}, tip[]}, 
               {Null, {Null, {Null, tip[], tip[], tip[]}, 
                 {Null, {Null, {Null, tip[], tip[]}, tip[]}, 
                  {Null, {Null, 
                    {Null, {Null, 
                      {Null, 
                       {Null, {Null, {Null, tip[], tip[]}, tip[]}, 
                       tip[]}, tip[]}, tip[]}, 
                     {Null, {Null, tip[], tip[]}, 
                      {Null, tip[], tip[], tip[]}}}, 
                    {Null, {Null, tip[], tip[]}, tip[]}}, 
                   {Null, tip[], tip[]}}, tip[]}}, 
                {Null, {Null, tip[], tip[]}, tip[]}}}}}}}}, tip[]}, 
        {Null, tip[], {Null, {Null, tip[], {Null, tip[], tip[]}}, 
          {Null, {Null, {Null, tip[], tip[]}, 
            {Null, {Null, tip[], tip[]}, 
             {Null, tip[], {Null, tip[], tip[]}}}}, 
           {Null, {Null, tip[], {Null, tip[], tip[]}}, 
            {Null, {Null, {Null, tip[], {Null, tip[], tip[]}}, 
              {Null, tip[], tip[]}}, 
             {Null, {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}, 
              {Null, tip[], {Null, {Null, tip[], tip[]}, 
                {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}}}}}}}}}}, 
       tip[], {Null, {Null, {Null, tip[], {Null, tip[], tip[]}}, 
         {Null, tip[], tip[]}}, 
        {Null, {Null, tip[], {Null, {Null, tip[], {Null, tip[], tip[]}}, 
           {Null, {Null, {Null, {Null, tip[], tip[]}, tip[]}, 
             {Null, tip[], tip[]}}, 
            {Null, {Null, tip[], tip[]}, 
             {Null, {Null, {Null, tip[], {Null, tip[], tip[]}}, 
               {Null, tip[], tip[]}, {Null, tip[], tip[], tip[], tip[]}}, 
              {Null, tip[], tip[]}}}}}}, 
         {Null, {Null, {Null, tip[], {Null, tip[], {Null, tip[], tip[]}}}, 
           {Null, tip[], tip[]}, {Null, tip[], tip[], tip[]}, 
           {Null, tip[], tip[]}}, {Null, tip[], tip[]}}}}, 
       {Null, tip[], {Null, tip[], 
         {Null, {Null, tip[], {Null, {Null, tip[], tip[]}, 
            {Null, tip[], {Null, {Null, tip[], tip[]}, 
              {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}, tip[]}}}}}, 
          {Null, tip[], {Null, 
            {Null, tip[], {Null, tip[], 
              {Null, tip[], {Null, {Null, tip[], tip[]}, 
                {Null, {Null, tip[], {Null, tip[], tip[]}}, 
                 {Null, tip[], {Null, tip[], {Null, tip[], tip[]}}}}}}}}, 
            {Null, {Null, {Null, {Null, tip[], tip[]}, tip[]}, tip[]}, 
             {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}}, 
            {Null, {Null, {Null, tip[], tip[]}, 
              {Null, {Null, tip[], tip[], {Null, tip[], tip[], tip[]}}, 
               {Null, tip[], {Null, tip[], 
                 {Null, {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}, 
                  {Null, tip[], tip[]}}}}}}, 
             {Null, {Null, tip[], {Null, tip[], tip[]}, 
               {Null, tip[], tip[]}}, 
              {Null, {Null, {Null, tip[], tip[]}, tip[]}, 
               {Null, tip[], {Null, tip[], tip[]}, 
                {Null, {Null, {Null, tip[], tip[]}, tip[]}, 
                 {Null, {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}, 
                  tip[]}}, {Null, 
                 {Null, tip[], 
                  {Null, {Null, {Null, tip[], tip[]}, tip[]}, tip[]}, 
                  {Null, {Null, 
                    {Null, tip[], 
                     {Null, tip[], {Null, tip[], tip[]}, 
                      {Null, tip[], tip[]}}}, tip[]}, tip[]}}, 
                 {Null, {Null, 
                   {Null, {Null, tip[], tip[]}, {Null, tip[], tip[]}}, 
                   tip[]}, {Null, tip[], 
                   {Null, {Null, {Null, {Null, tip[], tip[]}, tip[]}, 
                     {Null, {Null, {Null, {Null, tip[], tip[]}, tip[]}, 
                       {Null, tip[], tip[], tip[]}}, tip[]}}, 
                    {Null, {Null, tip[], tip[]}, 
                     {Null, {Null, {Null, {Null, tip[], tip[]}, tip[]}, 
                       tip[]}, tip[]}}}}}}}}}}}}}}}}}}}}};