#NEXUS
BEGIN TREES;
    tree tree_from_phyml = ((((((((Chimp:0.002646154700000014,Human:0.0025360128999999954):0.004974031200000006,Orangutan:0.007858024899999982):0.005083428599999995,Macaque:0.0244279013):0.03277411450000001,(Bushbaby:0.0584035061,Mouselemur:0.035090660900000015):0.012689592199999988):0.006351216899999956,((((Rat:0.04082554840000002,Mouse:0.03594933789999999):0.09737204719999998,Guinea_pig:0.0959334933):0.003419646600000048,G_squirrel:0.06196620829999999):0.008530266299999978,(Rabbit:0.0449045213,Pika:0.08344542650000003):0.037343581999999986):0.00420666550000004):0.007465057999999969,((Shrew:0.12675756719999998,Hedgehog:0.09829135979999998):0.018543300799999995,(Cat:0.03768168839999997,Dog:0.03889614089999999):0.02021927229999998):0.007764776000000029):0.008090884500000006,(Elephant:0.045532129899999996,Tenrec:0.1130345391):0.020299659999999997):0.05398220695,(Platypus:0.2590121974,Opossum:0.1674570976):0.05398220695);
END;

begin r8s;

  [ # HOW TO DO A PROPER ANALYSIS? ]
  [ # 1. get the best tree with branch lengths ]
  [ # 2. get the best calibration points ]
  [ # 3. CROSSVALIDATION for simple methods LF and NPRS with all algorithms ]
  [ # 4. CROSSVALIDATION for PL with different smoothing for all algoraithms and penalty log and add ]
  [ # 5. Select the method with lowest chi-squared error. ]
  [ #5.1 use directly TN which is the one the author recommends. ]
  [ # 6. Do the analysis with selected method, algorithm (and smoothing) and check_gradien ] 

 [ #set the number of sites since the tree is usually in per site units i.e. MrBayes or Phyml.]
 BLFORMAT nsites= 100649 lengths=persite ultrametric=no; 
 
 [# collapse;][ # collapse 0 length branches, polytomies are okay]
 
 [ #prune taxon= Shrew;][ # drop a taxa we don't want to use, this is usually important as you root with it in tree building but don't ]

[ # set some ages for internal nodes (which were labeled in the tree file, you can also define an internal label with the mrca command)]
 MRCA A Opossumhis hippopotamus;
 fixage taxon=A age=1; 

 [ # The fidelity with which any of the three methods explain the branch length variation can be explored using a cross-validation option ]
    set ftol=1e-9 maxiter=2000; [ #this line doesnt change results for TN]
    [DIVTIME crossv=yes method=LF algorithm=TN ;
    DIVTIME crossv=yes method=LF algorithm=POWELL;
    DIVTIME crossv=yes method=NPRS algorithm=POWELL;
    DIVTIME crossv=yes method=LF algorithm=QNEWT;]
   

[ #With PL we tried to optimize ‘‘smoothing’’ to correct for autocorrelation (R8S,
  #v1.70, user’s manual, available from M. Sanderson at http://ginger.ucdavis.edu/r8s/). The optimization involved
  #varying the value of the smoothing parameter, at increments of 0.5, from 22 to 9 on a log10 scale and, at each
  #value for the smoothing parameter, having the program estimate the prediction error that resulted from pruning
  #each terminal branch in turn. In our dataset, the prediction error did not show a minimum; instead it continued to
  #decrease even when the value of the smoothing parameter became high. This lack of a minimum indicated that
  #substitution rates were not correlated by lineage, but they could have been either clocklike or chaotic (R8S, v1.70,
  #user’s manual, available from M. Sanderson at http://ginger.ucdavis.edu/r8s/). Because the horizontal root-to-tip
  #branch lengths in the tree (FIG. 2) varied, we knew that the rates had not been clock-like. This led us to conclude that
  #rates had been chaotic. In the absence of predictable, lineage-specific rate correlations, we followed the recommendation
  #of the R8S program documentation and analyzed the data with LF, assuming a global substitution rate.
  #]
  [ #Repeat analysis Taylor and Barbee to determine what to use PL or LF]
  divtime penalty=log method=pl algorithm=tn crossv=yes cvstart=-8 cvinc=0.5 cvnum=8;
  [divtime penalty=add method=pl algorithm=tn crossv=yes cvstart=-4.65 cvinc=0.01 cvnum=10;]
  [divtime penalty=log method=pl algorithm=powell crossv=yes cvstart=-4.65 cvinc=0.01 cvnum=10;]
  [divtime penalty=add method=pl algorithm=powell crossv=yes cvstart=-4.65 cvinc=0.01 cvnum=10;]
  [divtime penalty=log method=pl algorithm=qnewt crossv=yes cvstart=-0.5 cvinc=0.1 cvnum=10;]
  [divtime penalty=add method=pl algorithm=qnewt crossv=yes cvstart=-0.5 cvinc=0.1 cvnum=10;]

  [ #Look for the smallest value in the last ("Chi Square Error") column, then choose the value from the second ("smooth") column in that same row.]
  [ # EXAMPLE ]
  [ #NEXUS ]
  [ #collapse; ]
  [ #fixage taxon=lp age=420; ]
  [ #set ftol=1e-7 maxiter=2000; ]
  [ #set verbose=1; [suppresses huge amount of output in CV analyses] ]
  [ #[divtime method=pl algorithm=tn cvStart=0 cvInc=0.5 cvNum=8 crossv=yes;] ]
  [ #set smoothing=2.5; ]
  [ #divtime method=pl algorithm=tn; ]
  [ #describe plot=tree_description; ]
  [ #end; ]
   

   [ #We tested for possible alternative
  # equally likely estimates for node ages with the LF program’s
  # CHECKGRADIENT option and set the analysis for 10 replicates
  # under different starting conditions, using ‘‘set num_time_-
  # guesses’’. 
  # The command set seed=k, where k is some integer, should be invoked to initialize the random number generator.
  ]

 [ #estimate divergence time under the langey-fitch method using a powell algorithm]
 [ set smoothing=2.5e-05;
  set checkGradient=yes num_time_guesses=10 seed=123456789;
  set penalty=log;
  DIVTIME method=PL algorithm=TN ;]
 

 [ #print out the table with ages for all the internal nodes]
 [showage;
 describe plot=chronogram;
 describe plot=chrono_description;
]


[ ###### How to run the file
  #Batch mode is invoked similarly but with the addition of the –b switch:
  # r8s –b –f datafile <cr>]


end;