1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
|
/* This is an example HY-PHY Batch File.
In this file, we will illustrate how to use the "ReplicateConstraint" command to
apply a constraint to the tree branches.
This example used "ChoiceList" to obtain user input, and "ReplicateConstraint"
to set up a simple LRT test
Sergei L. Kosakovsky Pond and Spencer V. Muse
June 2001.
*/
/* 1. Include a file which reads the data and the tree,
defines a model and the likelihood function */
#include "shared.bf";
/* 2. We begin by obtaining the MLE estimates for the
tree with each branch having 2 rates - transitions
and transversions, and printing out the result.
We will need this for the LRT later
*/
Optimize (res, hivLik);
fprintf (stdout, "1).Unconstrained HKY85\n\n", hivLik);
/* 3. Use ChoiceList to obtain user input on what model
parameters to constrain.
*/
ChoiceList (userChoice,"Parameter to Constrain",1,SKIP_NONE,
"Transitions","Constrain transition rates",
"Transversions","Constrain transversion rates",
"Both","Constraint transition and transversion rates");
/* 4. If the selection was canceled "userChoice" will be set to -1,
otherwise it will contain the index of the selection, starting at 0.
*/
if (userChoice<0)
{
return;
}
/* 5. Now we set up the constraint string to pass to ReplicateConstraint. */
if (userChoice==0)
{
constraintString = "trst";
}
else
{
if (userChoice==1)
{
constraintString = "trvs";
}
else
{
if (userChoice==2)
{
constraintString = "?";
}
}
}
constraintString = "this1.?."+constraintString+":=this2.?."+constraintString;
/* 6. Set up the constraint on subtrees starting at iNode1 and iNode2 and obtain the MLE */
ReplicateConstraint (constraintString,hivTree.iNode1,hivTree.iNode2);
Optimize (res2, hivLik);
fprintf (stdout, "\n\n2).Constrained model\n\n", hivLik);
/* 7. Perform the LRT */
lrStat = 2(res[1][0]-res2[1][0]);
degFreedom = (res[1][1]-res2[1][1]);
fprintf (stdout, "\n\n LRT:\n\n lrt statistic = ",lrStat,"\n degrees of freedom = ", degFreedom,
"\n\n Chi-squared P-Value = ", 1-CChi2(lrStat,degFreedom));
/* Note: take a peek at RelativeRate.bf in TemplateBatchFiles to
see how to use the ChoiceList in case the list of available model
parameters is not known apriori */
|
