Skip to content

Commit ed42400

Browse files
aarmeyclaude
aarmey
and
claude
committed
Add BDMM-Prime comparison inputs and results to the KP-Tracer report
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01A7EXGcBoM1ESZnrruqPQff
1 parent 6461b95 commit ed42400

7 files changed

Lines changed: 576 additions & 7 deletions

File tree

‎analysis/kptracer/REPORT.md‎

Lines changed: 18 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -80,8 +80,24 @@ Observations:
8080

8181
**Comparison with BDMM-Prime on a 100-cell tree** (`bdmm_comparison.py`, `results/bdmm_thmm.json`). Both methods get the same simulated 100-leaf time tree, the true tip types (4 types), and the same branch lengths.
8282

83-
- **tHMM:** fits in 0.4 s, with emissions fixed at near-indicators of the tip type. A 200-replicate parametric bootstrap takes 165 s, and its 95% intervals cover all 12 true rates. The four large rates are recovered in the right order: 0→1 at 0.69 (true 0.45), 1→0 at 0.27 (0.15), 1→2 at 0.26 (0.30), 2→3 at 0.14 (0.24). With only 34 true switches on the tree, the small rates go to the zero boundary.
84-
- **BDMM-Prime:** BDMM_RESULTS_PLACEHOLDER
83+
- **tHMM:** fits in 0.4 s, with emissions fixed at near-indicators of the tip type. A 200-replicate parametric bootstrap takes 165 s. Its 95% intervals cover all 12 true rates, but they are very wide on a tree this small. The four large rates are recovered in the right order: 0→1 at 0.69 (true 0.45), 1→0 at 0.27 (0.15), 1→2 at 0.26 (0.30), 2→3 at 0.14 (0.24). With only 34 true switches on the tree, the small rates go to the zero boundary.
84+
- **BDMM-Prime** (BEAST 2.7.7, BDMM-Prime 2.7.2; XML and summary in `bdmm/`):
85+
- **Setup:** fixed tree; the standard BDMM likelihood integrated over type histories; shared birth rate; no death; ρ = 1; Exponential(mean 0.3) priors on the 12 type-change rates; estimated root frequencies.
86+
- **Runtime:** about 2.2–2.6 ms per MCMC state. Mixing was excellent (single-chain ESS > 2,400 over 2.5M states), so ESS > 200 needs about 200k states, roughly 7–10 minutes for one chain. The full 4-chain run took 107 minutes.
87+
- **Accuracy:** its 95% HPDs cover all 12 true rates, but they are wide. The rates out of the rare type 3 sit near the prior median (0.21–0.27 against true 0.015–0.03), and the root type is not identified.
88+
89+
| rate | true | tHMM MLE [bootstrap 95%] | BDMM-Prime median [95% HPD] |
90+
|---|---|---|---|
91+
| 0→1 | 0.45 | 0.69 [0.01, 3.1] | 0.51 [0.10, 0.95] |
92+
| 1→0 | 0.15 | 0.27 [0.00, 2.5] | 0.29 [0.04, 0.59] |
93+
| 1→2 | 0.30 | 0.26 [0.00, 3.4] | 0.20 [0.03, 0.41] |
94+
| 2→3 | 0.24 | 0.14 [0.00, 2.7] | 0.13 [0.00, 0.38] |
95+
| rates out of type 3 | 0.015–0.03 | 0.00–0.01 | 0.21–0.27 (prior-driven) |
96+
97+
On the same information the two methods are about equally accurate for the well-determined rates. The tHMM is about 1,000× faster than even a minimal BDMM-Prime run. It also scales linearly to the 14k-cell KP-Tracer trees, where BDMM-type inference is not practical. The trade-offs:
98+
- The tHMM conditions on the tree and does not model type-dependent birth or death, which BDMM does.
99+
- Its intervals come from a parametric bootstrap, not a posterior. On 100 cells they are much wider than the HPDs, because each refit of a 12-rate model to about 34 switches occasionally explodes one rate, and nothing like BDMM-Prime's prior holds rates down.
100+
- Its maximum-likelihood rates for rarely visited states sit on the zero boundary, where BDMM-Prime's are prior-driven.
85101

86102
## KP-Tracer analysis
87103

‎analysis/kptracer/bdmm/bdmm_fixedtree.xml‎

Lines changed: 281 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,281 @@
1+
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
2+
<beast version="2.7" namespace="beast.base.inference:beast.base.inference.parameter:beast.base.inference.distribution:beast.base.inference.operator:beast.base.evolution.tree:beast.base.evolution.alignment:beast.base.evolution.operator:beast.base.inference.util:bdmmprime.distribution:bdmmprime.parameterization:bdmmprime.util">
3+
4+
<taxonset id="taxa" spec="TaxonSet">
5+
<taxon id="n123" spec="Taxon"/>
6+
<taxon id="n185" spec="Taxon"/>
7+
<taxon id="n186" spec="Taxon"/>
8+
<taxon id="n130" spec="Taxon"/>
9+
<taxon id="n66" spec="Taxon"/>
10+
<taxon id="n175" spec="Taxon"/>
11+
<taxon id="n176" spec="Taxon"/>
12+
<taxon id="n159" spec="Taxon"/>
13+
<taxon id="n160" spec="Taxon"/>
14+
<taxon id="n137" spec="Taxon"/>
15+
<taxon id="n33" spec="Taxon"/>
16+
<taxon id="n34" spec="Taxon"/>
17+
<taxon id="n149" spec="Taxon"/>
18+
<taxon id="n150" spec="Taxon"/>
19+
<taxon id="n47" spec="Taxon"/>
20+
<taxon id="n96" spec="Taxon"/>
21+
<taxon id="n100" spec="Taxon"/>
22+
<taxon id="n99" spec="Taxon"/>
23+
<taxon id="n197" spec="Taxon"/>
24+
<taxon id="n198" spec="Taxon"/>
25+
<taxon id="n177" spec="Taxon"/>
26+
<taxon id="n178" spec="Taxon"/>
27+
<taxon id="n191" spec="Taxon"/>
28+
<taxon id="n192" spec="Taxon"/>
29+
<taxon id="n139" spec="Taxon"/>
30+
<taxon id="n181" spec="Taxon"/>
31+
<taxon id="n182" spec="Taxon"/>
32+
<taxon id="n131" spec="Taxon"/>
33+
<taxon id="n106" spec="Taxon"/>
34+
<taxon id="n127" spec="Taxon"/>
35+
<taxon id="n121" spec="Taxon"/>
36+
<taxon id="n122" spec="Taxon"/>
37+
<taxon id="n135" spec="Taxon"/>
38+
<taxon id="n136" spec="Taxon"/>
39+
<taxon id="n143" spec="Taxon"/>
40+
<taxon id="n144" spec="Taxon"/>
41+
<taxon id="n195" spec="Taxon"/>
42+
<taxon id="n196" spec="Taxon"/>
43+
<taxon id="n147" spec="Taxon"/>
44+
<taxon id="n148" spec="Taxon"/>
45+
<taxon id="n45" spec="Taxon"/>
46+
<taxon id="n46" spec="Taxon"/>
47+
<taxon id="n37" spec="Taxon"/>
48+
<taxon id="n165" spec="Taxon"/>
49+
<taxon id="n110" spec="Taxon"/>
50+
<taxon id="n171" spec="Taxon"/>
51+
<taxon id="n172" spec="Taxon"/>
52+
<taxon id="n69" spec="Taxon"/>
53+
<taxon id="n71" spec="Taxon"/>
54+
<taxon id="n193" spec="Taxon"/>
55+
<taxon id="n194" spec="Taxon"/>
56+
<taxon id="n187" spec="Taxon"/>
57+
<taxon id="n188" spec="Taxon"/>
58+
<taxon id="n102" spec="Taxon"/>
59+
<taxon id="n103" spec="Taxon"/>
60+
<taxon id="n104" spec="Taxon"/>
61+
<taxon id="n53" spec="Taxon"/>
62+
<taxon id="n85" spec="Taxon"/>
63+
<taxon id="n86" spec="Taxon"/>
64+
<taxon id="n88" spec="Taxon"/>
65+
<taxon id="n125" spec="Taxon"/>
66+
<taxon id="n73" spec="Taxon"/>
67+
<taxon id="n151" spec="Taxon"/>
68+
<taxon id="n152" spec="Taxon"/>
69+
<taxon id="n77" spec="Taxon"/>
70+
<taxon id="n78" spec="Taxon"/>
71+
<taxon id="n120" spec="Taxon"/>
72+
<taxon id="n133" spec="Taxon"/>
73+
<taxon id="n134" spec="Taxon"/>
74+
<taxon id="n141" spec="Taxon"/>
75+
<taxon id="n142" spec="Taxon"/>
76+
<taxon id="n183" spec="Taxon"/>
77+
<taxon id="n184" spec="Taxon"/>
78+
<taxon id="n117" spec="Taxon"/>
79+
<taxon id="n118" spec="Taxon"/>
80+
<taxon id="n60" spec="Taxon"/>
81+
<taxon id="n91" spec="Taxon"/>
82+
<taxon id="n92" spec="Taxon"/>
83+
<taxon id="n161" spec="Taxon"/>
84+
<taxon id="n157" spec="Taxon"/>
85+
<taxon id="n158" spec="Taxon"/>
86+
<taxon id="n167" spec="Taxon"/>
87+
<taxon id="n168" spec="Taxon"/>
88+
<taxon id="n63" spec="Taxon"/>
89+
<taxon id="n113" spec="Taxon"/>
90+
<taxon id="n179" spec="Taxon"/>
91+
<taxon id="n180" spec="Taxon"/>
92+
<taxon id="n108" spec="Taxon"/>
93+
<taxon id="n163" spec="Taxon"/>
94+
<taxon id="n164" spec="Taxon"/>
95+
<taxon id="n84" spec="Taxon"/>
96+
<taxon id="n89" spec="Taxon"/>
97+
<taxon id="n90" spec="Taxon"/>
98+
<taxon id="n153" spec="Taxon"/>
99+
<taxon id="n154" spec="Taxon"/>
100+
<taxon id="n155" spec="Taxon"/>
101+
<taxon id="n156" spec="Taxon"/>
102+
<taxon id="n170" spec="Taxon"/>
103+
<taxon id="n189" spec="Taxon"/>
104+
<taxon id="n190" spec="Taxon"/>
105+
</taxonset>
106+
107+
<!-- FIXED tree: not in the state, no tree operators -->
108+
<tree id="tree" spec="TreeParser" IsLabelledNewick="true" adjustTipHeights="true" taxonset="@taxa"
109+
newick="((((n185:0.105435,n186:0.105435)n75:0.758506,((n175:0.169748,n176:0.169748)n145:0.112940,(n159:0.234216,n160:0.234216)n146:0.048472)n76:0.581254)n43:0.449058,(n123:0.450563,((n137:0.324692,(n149:0.272999,n150:0.272999)n138:0.051694)n129:0.099779,n130:0.424471)n124:0.026092)n44:0.862437)n1:1.968832,(((((n47:1.159471,(n181:0.117813,n182:0.117813)n48:1.041658)n23:1.098771,((n131:0.409898,(n165:0.217494,(n193:0.047186,n194:0.047186)n166:0.170309)n132:0.192403)n55:0.683481,(((n187:0.104629,n188:0.104629)n109:0.501520,n110:0.606149)n105:0.017891,n106:0.624040)n56:0.469339)n24:1.164863)n15:0.438788,(((n127:0.431499,(n171:0.173877,n172:0.173877)n128:0.257622)n95:0.332616,n96:0.764115)n79:0.093977,(n100:0.734049,n99:0.734049)n80:0.124044)n16:1.838938)n5:0.300383,((n33:1.783769,n34:1.783769)n21:0.587406,(((n121:0.490003,n122:0.490003)n111:0.102243,(n135:0.369630,n136:0.369630)n112:0.222616)n97:0.169138,(n197:0.012586,n198:0.012586)n98:0.748797)n22:1.609791)n6:0.626239)n3:0.015266,((((n177:0.135231,n178:0.135231)n173:0.037077,(n191:0.052557,n192:0.052557)n174:0.119751)n65:0.745767,n66:0.918074)n7:2.050193,(((n139:0.323601,(n143:0.288174,n144:0.288174)n140:0.035427)n41:1.038820,((n195:0.015571,n196:0.015571)n115:0.514336,(n147:0.281240,n148:0.281240)n116:0.248668)n42:0.832513)n10:1.591565,((((n69:0.909870,((n125:0.436570,(n183:0.108761,n184:0.108761)n126:0.327809)n101:0.212951,n102:0.649521)n70:0.260349)n19:1.594987,((n103:0.645152,n104:0.645152)n49:0.464789,(n53:1.095520,(n73:0.887923,((n157:0.235615,n158:0.235615)n81:0.607617,(n167:0.214877,n168:0.214877)n82:0.628355)n74:0.044692)n54:0.207596)n50:0.014421)n20:1.394916)n17:0.101622,((n71:0.903043,(n85:0.829271,n86:0.829271)n72:0.073772)n25:1.128496,(((n151:0.267928,n152:0.267928)n31:1.569555,(n77:0.859500,n78:0.859500)n32:0.977983)n27:0.069964,(((n117:0.510409,n118:0.510409)n35:1.142729,((n63:0.927585,((n153:0.264503,n154:0.264503)n107:0.346663,n108:0.611166)n64:0.316418)n59:0.114692,n60:1.042276)n36:0.610861)n29:0.219068,((n91:0.777811,n92:0.777811)n57:0.302458,((n113:0.531419,(n163:0.218871,n164:0.218871)n114:0.312548)n61:0.413352,((((n155:0.236850,n156:0.236850)n93:0.527531,((n189:0.086787,n190:0.086787)n169:0.097964,n170:0.184751)n94:0.579630)n83:0.066969,n84:0.831350)n67:0.078867,(n89:0.779814,n90:0.779814)n68:0.130402)n62:0.034554)n58:0.135498)n30:0.791937)n28:0.035241)n26:0.124093)n18:0.574939)n11:0.316059,((n45:1.179748,n46:1.179748)n13:1.605255,(n37:1.583666,((((n161:0.221276,(n179:0.121742,n180:0.121742)n162:0.099534)n119:0.287553,n120:0.508829)n51:0.592224,(n133:0.382097,n134:0.382097)n52:0.718956)n39:0.284872,((n141:0.313129,n142:0.313129)n87:0.509341,n88:0.822470)n40:0.563455)n38:0.197741)n14:1.201337)n12:0.137534)n9:0.031449)n8:0.014281)n4:0.044412)n2:0.269154)n0:0.000000;"/>
110+
111+
<typeTraitSet id="typeTraitSet" spec="TraitSet" traitname="type" taxa="@taxa">
112+
n123=1,
113+
n185=2,
114+
n186=1,
115+
n130=0,
116+
n66=1,
117+
n175=0,
118+
n176=0,
119+
n159=0,
120+
n160=0,
121+
n137=0,
122+
n33=0,
123+
n34=1,
124+
n149=0,
125+
n150=0,
126+
n47=0,
127+
n96=2,
128+
n100=2,
129+
n99=2,
130+
n197=2,
131+
n198=2,
132+
n177=0,
133+
n178=1,
134+
n191=1,
135+
n192=1,
136+
n139=2,
137+
n181=0,
138+
n182=0,
139+
n131=2,
140+
n106=1,
141+
n127=2,
142+
n121=2,
143+
n122=2,
144+
n135=2,
145+
n136=2,
146+
n143=2,
147+
n144=2,
148+
n195=1,
149+
n196=1,
150+
n147=1,
151+
n148=1,
152+
n45=2,
153+
n46=1,
154+
n37=2,
155+
n165=2,
156+
n110=1,
157+
n171=3,
158+
n172=3,
159+
n69=1,
160+
n71=1,
161+
n193=2,
162+
n194=2,
163+
n187=1,
164+
n188=1,
165+
n102=0,
166+
n103=1,
167+
n104=1,
168+
n53=2,
169+
n85=1,
170+
n86=1,
171+
n88=0,
172+
n125=0,
173+
n73=0,
174+
n151=3,
175+
n152=3,
176+
n77=2,
177+
n78=1,
178+
n120=1,
179+
n133=0,
180+
n134=1,
181+
n141=1,
182+
n142=1,
183+
n183=0,
184+
n184=0,
185+
n117=1,
186+
n118=1,
187+
n60=1,
188+
n91=1,
189+
n92=0,
190+
n161=1,
191+
n157=0,
192+
n158=0,
193+
n167=0,
194+
n168=0,
195+
n63=2,
196+
n113=0,
197+
n179=1,
198+
n180=1,
199+
n108=1,
200+
n163=0,
201+
n164=0,
202+
n84=0,
203+
n89=1,
204+
n90=1,
205+
n153=1,
206+
n154=1,
207+
n155=0,
208+
n156=1,
209+
n170=0,
210+
n189=1,
211+
n190=1
212+
</typeTraitSet>
213+
214+
<typeSet id="typeSet" spec="TypeSet" typeTraitSet="@typeTraitSet"/>
215+
<processLength id="processLength" spec="ProcessLength" tree="@tree"/>
216+
217+
<run id="mcmc" spec="MCMC" chainLength="2500000">
218+
<state id="state" storeEvery="5000">
219+
<stateNode id="birthRate" spec="RealParameter" lower="0.0" dimension="1">1.0</stateNode>
220+
<stateNode id="migrationRate" spec="RealParameter" lower="0.0" dimension="12">0.1</stateNode>
221+
<stateNode id="rootFreqs" spec="RealParameter" lower="0.0" upper="1.0" dimension="4">0.25</stateNode>
222+
</state>
223+
224+
<distribution id="posterior" spec="CompoundDistribution">
225+
<distribution id="prior" spec="CompoundDistribution">
226+
<distribution id="birthRatePrior" spec="Prior" x="@birthRate">
227+
<distr spec="LogNormalDistributionModel" meanInRealSpace="false" M="0.0" S="1.0"/>
228+
</distribution>
229+
<distribution id="migrationRatePrior" spec="Prior" x="@migrationRate">
230+
<distr spec="Exponential" mean="0.3"/>
231+
</distribution>
232+
<distribution id="rootFreqsPrior" spec="Prior" x="@rootFreqs">
233+
<distr spec="Dirichlet"><parameter spec="RealParameter" name="alpha" value="1.0 1.0 1.0 1.0"/></distr>
234+
</distribution>
235+
</distribution>
236+
<distribution id="likelihood" spec="CompoundDistribution">
237+
<distribution id="treePrior" spec="BirthDeathMigrationDistribution" tree="@tree"
238+
typeTraitSet="@typeTraitSet" startTypePriorProbs="@rootFreqs"
239+
conditionOnRoot="true" parallelize="true">
240+
<parameterization id="param" spec="CanonicalParameterization" typeSet="@typeSet" processLength="@processLength">
241+
<birthRate spec="SkylineVectorParameter" typeSet="@typeSet" processLength="@processLength" skylineValues="@birthRate"/>
242+
<deathRate spec="SkylineVectorParameter" typeSet="@typeSet" processLength="@processLength">
243+
<skylineValues spec="RealParameter" value="0.0"/>
244+
</deathRate>
245+
<samplingRate spec="SkylineVectorParameter" typeSet="@typeSet" processLength="@processLength">
246+
<skylineValues spec="RealParameter" value="0.0"/>
247+
</samplingRate>
248+
<removalProb spec="SkylineVectorParameter" typeSet="@typeSet" processLength="@processLength">
249+
<skylineValues spec="RealParameter" value="1.0"/>
250+
</removalProb>
251+
<rhoSampling spec="TimedParameter" typeSet="@typeSet" processLength="@processLength" timesAreAges="true">
252+
<times spec="RealParameter" value="0.0"/>
253+
<values spec="RealParameter" value="1.0"/>
254+
</rhoSampling>
255+
<migrationRate spec="SkylineMatrixParameter" typeSet="@typeSet" processLength="@processLength" skylineValues="@migrationRate"/>
256+
</parameterization>
257+
</distribution>
258+
</distribution>
259+
</distribution>
260+
261+
<operator id="birthRateScale" spec="ScaleOperator" parameter="@birthRate" scaleFactor="0.5" weight="3.0"/>
262+
<operator id="migRateScale" spec="ScaleOperator" parameter="@migrationRate" scaleFactor="0.5" weight="12.0"/>
263+
<operator id="migRateScaleAll" spec="ScaleOperator" parameter="@migrationRate" scaleAll="true" scaleFactor="0.75" weight="2.0"/>
264+
<operator id="rootFreqsExchange" spec="DeltaExchangeOperator" parameter="@rootFreqs" delta="0.2" weight="2.0"/>
265+
266+
<logger id="tracelog" spec="Logger" fileName="bdmm_fixedtree_chain1.log" logEvery="500" sanitiseHeaders="true">
267+
<log idref="posterior"/>
268+
<log idref="likelihood"/>
269+
<log idref="prior"/>
270+
<log idref="treePrior"/>
271+
<log idref="birthRate"/>
272+
<log idref="migrationRate"/>
273+
<log idref="rootFreqs"/>
274+
</logger>
275+
<logger id="screenlog" spec="Logger" logEvery="25000">
276+
<log idref="posterior"/>
277+
<log idref="likelihood"/>
278+
<log spec="ESS" arg="@posterior"/>
279+
</logger>
280+
</run>
281+
</beast>

0 commit comments

Comments
 (0)