Meta tags:
Headings (most frequently used words):
model, 2026, character, multi, july, discrete, dependent, ou, phytools, for, regime, in, about, blog, joint, under, the, using, on, phylogenetic, tools, comparative, biology, pages, sunday, august, wednesday, 29, tuesday, 14, monday, 13, saturday, 11, this, me, links, archive, small, experiment, but, strong, case, estimation, of, continuous, trait, evolution, stochastic, mapping, threshold, fitthresh, comparing, to, ouwie, visualizing, process, tree, more, optimum, null, ornstein, uhlenbeck,
Text of the page (most frequently used words):
the (244), model (132), this (93), our (92), discrete (91), theta (83), #character (82), alpha (70), for (65), that (65), and (61), with (55), low (54), high (53), phytools (49), tree (48), phy (47), trait (46), log (45), set (45), letters (42), sigsq (41), using (39), history (38), let (37), multi (37), levels (36), ouwie (34), fit (33), setnames (32), but (31), fitmultiou (31), continuous (31), true (30), sim_tree (30), likelihood (29), 100 (29), regime (29), fitted (28), data (28), threshold (27), 200 (26), now (25), object (24), fitmk (24), root (23), simulate (23), plot (23), matrix (22), colors (22), smp (22), stochastic (21), july (20), going (20), just (20), may (19), some (19), can (19), share (18), mk_fit (18), parameters (18), tips (18), lwd (18), not (18), 267 (18), phylogenetic (17), parameter (17), estimates (17), generating (17), hcl (17), includes (17), process (17), june (16), joint (16), levs (16), values (16), done (16), null (16), same (16), which (16), value (16), head (16), from (16), will (16), cols (16), december (15), january (15), march (15), august (15), september (15), fit_null (15), sig2 (15), has (15), solution (15), t49 (15), off (15), nodes (15), tip (15), dependent (15), init (15), april (14), november (14), more (14), class (14), ncores (14), sigma (14), simmap (14), october (13), should (13), estimated (13), prior (13), optimization (13), fitcontinuous (13), state (13), under (13), each (13), would (13), have (13), very (13), all (13), see (13), february (12), new (12), pretty (12), these (12), thinks (12), warning (12), geiger (12), multiou (12), simulation (12), use (12), ftype (12), internal (12), follows (12), max (12), merged (12), zero (12), t38 (12), rates (11), based (11), found (11), mle (11), get (11), states (11), upwards (11), col (11), labels (11), rooted (11), branch (11), length (11), because (11), fit_mou (11), fitthresh (11), t16 (11), t21 (11), t123 (11), t124 (11), loglik (10), compare (10), ou_fit (10), though (10), t19 (10), need (10), next (10), two (10), direction (10), also (10), transition (10), library (10), out (10), mean (10), map (10), results (10), ouwie_results (10), smps (10), revell (9), package (9), posted (9), here (9), optimizing (9), was (9), summary (9), iterations (9), function (9), start (9), t53 (9), pull (9), getstates (9), different (9), are (9), make (9), fixed (9), rstate (9), sum (9), probabilities (9), slightly (9), setting (9), other (9), liam (8), about (8), email (8), cool (8), null_logl (8), again (8), list (8), iter (8), parallel (8), trace (8), then (8), method (8), only (8), rep (8), par (8), legend (8), bottomleft (8), bty (8), mapped (8), load (8), nodeheights (8), across (8), maps (8), rate (8), models (8), case (7), 2026 (7), blog (7), expect (7), attr (7), even (7), closely (7), discretization (7), post (7), used (7), aic (7), number (7), might (7), sim (7), approximation (7), 283 (7), optimum (7), between (7), regimes (7), numeric (7), ms_thresh (7), liability (7), foreach (7), pinterest (6), facebook (6), blogthis (6), comments (6), 300 (6), comparative (6), best (6), factor (6), step (6), conditions (6), since (6), cex (6), lengths (6), anc (6), dimnames (6), pbtree (6), scale (6), time (6), how (6), close (6), 2094242 (6), 500 (6), 255589 (6), what (6), 2554668 (6), when (6), generate (6), changes (6), quite (6), genus_species (6), reg (6), t35 (6), frame (6), diag (6), t28 (6), t31 (6), t32 (6), t23 (6), t24 (6), mergemappedstates (6), liabilities (6), following (6), 472 (6), ouwie_data (6), small (5), posts (5), similar (5), 1883095 (5), well (5), equal (5), due (5), treating (5), nlminb (5), aicc (5), search (5), add (5), lend (5), estimation (5), trees (5), been (5), evolution (5), single (5), one (5), diffusion (5), print (5), 287 (5), seems (5), 400 (5), maxit (5), 2000 (5), run (5), first (5), already (5), false (5), course (5), fitting (5), t89 (5), t90 (5), vary (5), mapping (5), thresh_fit (5), sapply (5), 350 (5), histories (5), 2016 (4), you (4), before (4), version (4), both (4), likelihoods (4), increase (4), far (4), finally (4), once (4), opt (4), comparison (4), match (4), those (4), null_model (4), while (4), 188436 (4), free (4), t10 (4), t54 (4), simulated (4), specify (4), good (4), nrow (4), probably (4), note (4), them (4), few (4), pre (4), there (4), 7812 (4), better (4), than (4), 0693 (4), 600 (4), var (4), bit (4), starting (4), t44 (4), t57 (4), t136 (4), t177 (4), analysis (4), into (4), loading (4), thought (4), hidden (4), characters (4), original (4), infrequently (4), oum (4), finished (4), station (4), 02713072 (4), 142 (4), long (4), 028064 (4), 02921177 (4), they (4), finite (4), space (4), mkn (4), fitzjohn (4), nsim (4), mk_model (4), mkm (4), fun (4), bias (4), 7823842 (4), options (4), scipen (4), 355 (4), foo (4), ouwie_fits (4), type (4), 458785 (4), 250 (4), 03312624 (4), 241 (4), 036172 (4), 2013 (3), 2014 (3), experiment (3), page (3), tools (3), analyses (3), recommend (3), cran (3), home (3), 5433 (3), converge (3), obj (3), don (3), ready (3), takes (3), earlier (3), flat (3), convergence (3), diagnostics (3), failed (3), frequency (3), lik (3), bnd (3), bounds (3), res (3), iteration (3), maximum (3), still (3), root_state (3), code (3), actually (3), ahead (3), decided (3), reasonable (3), obtain (3), estimate (3), another (3), ornstein (3), uhlenbeck (3), boucher (3), démery (3), interested (3), biorxiv (3), weight (3), 0493 (3), know (3), much (3), 7939 (3), 0513 (3), 0689 (3), 4188483 (3), 700 (3), 06931472 (3), sensible (3), 284 (3), obtained (3), 511178 (3), 333333 (3), confirm (3), 0829 (3), 5109335 (3), 199 (3), want (3), work (3), recent (3), recently (3), readers (3), line (3), stat_theta (3), seq (3), min (3), packages (3), usually (3), way (3), begin (3), known (3), indeed (3), 8523 (3), 3246 (3), 4163 (3), 0271 (3), matches (3), ou1 (3), result (3), names (3), goes (3), great (3), rowsums (3), jointly (3), along (3), something (3), methods (3), biased (3), 0253 (3), 1673 (3), 352 (3), 354 (3), likely (3), 249 (3), 9258 (3), 159 (2), 2017 (2), 2025 (2), strong (2), lab (2), webpage (2), latest (2), similarity (2), 58921 (2), ignore (2), unlist (2), colnames (2), quick (2), too (2), 5912 (2), 5499 (2), 1496 (2), 430 (2), 8024 (2), evolutionary (2), traits (2), however (2), including (2), who (2), simulating (2), isn (2), demo (2), practice (2), packageversion (2), check (2), corresponding (2), trying (2), dependence (2), hope (2), separately (2), identifiable (2), without (2), ago (2), details (2), general (2), compared (2), anova (2), naturally (2), whether (2), were (2), 1335 (2), 1852 (2), 0707 (2), 2202 (2), 800 (2), 2003 (2), 900 (2), 2156 (2), 7829 (2), 2094 (2), converged (2), try (2), 1910 (2), 0528 (2), 2053 (2), 0527 (2), 9134 (2), fact (2), caution (2), 185 (2), 5743 (2), 0141 (2), unknown (2), taxa (2), none (2), hand (2), required (2), consider (2), future (2), morning (2), layout (2), mar (2), xlim (2), mtext (2), adj (2), transparent (2), predicted (2), distribution (2), stat (2), norm (2), singleton (2), later (2), illustration (2), look (2), fairly (2), keeping (2), mind (2), lnl (2), bic (2), ntax (2), 3007500 (2), 1198701 (2), optima (2), half (2), life (2), reporting (2), 304729 (2), arrived (2), reliable (2), fitoum (2), initializing (2), thorough (2), summarizing (2), algorithm (2), specified (2), defaulting (2), computing (2), determinant (2), inversion (2), vcv (2), 4519 (2), 1141 (2), 4520 (2), matter (2), got (2), fitou (2), give (2), simply (2), 3416784 (2), 6347880 (2), 3128115 (2), 5001694 (2), 7747887 (2), 5183583 (2), special (2), format (2), higher (2), fix (2), 6338 (2), 1838 (2), 175 (2), towards (2), infty (2), squared (2), sample (2), over (2), thesis (2), interesting (2), uses (2), review (2), comparing (2), thresholds (2), does (2), lowermost (2), through (2), complicated (2), many (2), back (2), makes (2), sense (2), standard (2), clearly (2), side (2), big (2), discretized (2), totally (2), times (2), t25 (2), t26 (2), t27 (2), t33 (2), thresh (2), underlying (2), functions (2), making (2), mvmorph (2), lot (2), days (2), away (2), seem (2), traditional (2), perhaps (2), difference (2), underestimate (2), 4625 (2), 4682 (2), 2082 (2), 473 (2), colmeans (2), average (2), 360 (2), 356 (2), 358 (2), printing (2), round (2), digits (2), rownames (2), makecluster (2), psock (2), registerdoparallel (2), optimize (2), dopar (2), stopcluster (2), parallelize (2), 142639 (2), 634035 (2), 209871 (2), 537554 (2), 871243 (2), 072081 (2), specific (2), mfrow (2), nulo (2), nuisance (2), differ (2), mathbf (2), ylim (2), 1009 (2), 9043 (2), 2719 (2), 0331 (2), 227 (2), 193 (2), cluster (2), doparallel (2), 6548297 (2), 7918578 (2), 1488759 (2), 4502061 (2), 3959087 (2), 2302126 (2), sunday (2), workflow (2), biology (2), powered, blogger, 2010, 162, 2011, 2012, 123, 2015, 132, 154, 2018, 2019, 2020, 2021, 2022, 2023, 2024, archive, umass, boston, faculty, phylogenetics, links, view, complete, profile, web, chronicles, development, install, attempting, replicate, any, site, beta, release, subscribe, atom, older, soon, come, micro, experimentation, confirms, limit, 5498946, 5498956, 1521741, 1495582, 5928578, 5912444, 1884358, cbind, throw, 0991, 0802, 0040, 0662, 5097, 6013, 3368, 0960, 1504, 6002, 5436, 1250, 2330, 1476, 5919, 5487, 1494, 1888, 5434, 331, 1883, having, effect, mode, had, realized, matched, ran, 29465, 549895, 152174, 592858, 294564, 589128, 017699, quickly, 5425203, 2690402, 4526387, 2366645, 1525741, 9399830, people, like, adapt, necessary, anyway, visualize, could, take, larger, started, colleague, prototype, consequently, write, today, demonstrating, illustrate, extent, numerical, precision, update, inadvertently, x_0, geq, github, weeks, technical, approach, saturday, preliminary, tells, basically, evidence, falls, 582, 0985, 166496e, 547, 5624, 000000e, explains, independent, easy, remarkably, remember, definitely, 2324, 6486, 9372, 0436, 2556, 315, 4681, 5480, 6211, 3112, 1711, 0816, 1765, 270, 5138, 9475, 9025, 0608, 1907, 0725, 2221, 268, 0952, 8971, 7990, 0222, 1960, 2176, 9648, 8941, 8285, 0684, 1935, 0701, 9573, 8831, 8318, 0694, 1936, 2166, 9496, 8828, 8724, 0711, 1903, 0708, 2171, 9394, 8728, 8737, 0625, 1894, 0710, 2152, 9342, 8513, 8355, 9719, 0724, 2116, 9037, 8589, 8295, 9641, 0720, 2139, 9002, 1000, 8613, 8330, 9623, 0719, 2149, 8995, 1100, 8608, 8220, 9628, 0713, 2154, 8972, 1200, 8600, 8276, 9672, 2008, 0712, 8957, 1300, 8269, 0429, 0894, 1860, 0700, 2134, 8096, 1400, 8193, 0484, 1319, 1836, 0691, 2161, 7998, 1500, 8143, 0503, 1171, 1844, 0696, 2132, 7953, 1600, 7951, 0509, 1326, 1842, 0690, 2107, 1700, 1340, 1851, 2095, 1723, 23239962, 64856597, 93721136, 04357107, 25558924, difficult, suspicious, 3575, 3680, 8215, 2403, 0115, 2066, 309, 6850, 3023, 0148, 6594, 2779, 0615, 1947, 8126, 2191, 1288, 9057, 2837, 0611, 1861, 2140, 1372, 9328, 2794, 0554, 1914, 0423, 2355, 1808, 9946, 2713, 0499, 1976, 9893, 2229, 1901, 0192, 2716, 0514, 2043, 9538, 2210, 9799, 2748, 0519, 2054, 9183, 2172, 2083, 9723, 2747, 2056, 9137, 2164, 2088, 9738, 2751, 9135, 9724, 2753, 2052, 913, 2157, 2097, 9725, 2752, farther, apart, comfortable, undoubtedly, substantially, 285, 5048, 195, 760975, 012919, 081833, 655514, 743851, 487702, 610151, hard, 0601, 2682, 0177, 409, 2704, 4702, 0182, 0846, 2622, 2005, 5735, 0150, 2569, 0941, 279, 2555, awesome, vector, 9418267, 5788125, 2039455, 7955580, 6064443, 9103976, nearly, allows, multiple, assumes, constant, continous, empirical, shall, aside, move, forward, readily, ape, writing, versions, checking, utmost, nonetheless, works, simpler, alternative, covered, monday, comment, nothing, heights, margins, second, subplot, phenogram, axis, las, equilibrium, dnorm, sqrt, polygon, border, unbranching, era, limits, 1001, chunk, objects, become, useful, contrast, three, phytools_liam, pic, twitter, com, jrmalmqpbh, https, lplhskjoo6, tweet, included, derived, visualizing, tuesday, folks, hold, whenever, unambiguous, diminish, increases, implements, 47594, 95187, 37292, 107, 3726, 75307314, 5184579, 09854245, 1132707, forget, almost, never, genuinely, knew, always, 2784, 6184, 1283, 0281, 1417, 4558, 8575, 3260, 1157, 0268, 4242, 8537, 3254, 1143, 0267, 4172, 4531, 8535, 3252, 1145, 0269, 4170, 1142, 404, 27841916, 61842864, 12828760, 02806402, identify, blogging, remind, reader, hypothesizing, hopefully, 25048, 501, 751, 150, 3165, 524044e, 510672e, 6250224, 3348563, 45480781, why, compile, put, bring, picture, updated, closer, currently, implemented, really, 18378, 41855, compute, assure, ourselves, converging, 168076, 000000, 085106, 625022, 250476, 500953, 750953, 086, 9456, 0721, 1329, 0914, 6004, 6244, 0116, 0897, 3007, 6123, 0053, 0892, 0283, 2222, 271, 0004, 0860, 0292, degree, exactly, assumptions, twice, completeness, point, generator, according, dependentn, hypothesis, explore, integrates, uncertainty, highly, implement, popular, brian, meara, jeremy, beaulieu, described, wow, beautiful, repeat, estimable, otherwise, non, normally, centers, whole, instead, unitless, scaleless, doesn, relative, positions, annoying, interpret, lower, 7568, 341537, 847661, sequence, threshstate, reuse, things, little, instance, crossing, involves, switches, forth, perfect, capturing, form, seeing, looks, running, turn, merge, contains, finely, effectively, typically, accomplished, modifying, iterating, subsequent, steps, message, elements, safely, believe, element, feel, undertake, implicit, within, 04421, 08842, 1050167, 90152, 80305, 8949833, symmetric, distinguishable, 901524, proceed, t37, t79, t80, t29, t30, t46, t47, else, 9810162, 1667724, 9394853, 5019343, 2609373, 1475927, 4670823, 2442964, 7697432, 0546414, fastbm, simplest, example, call, yet, software, straightforward, automated, show, inasmuch, crossed, fractal, moves, didn, trust, technically, enable, sorts, actual, weird, question, regarding, your, bins, test, multivariate, shape, reconstructing, ancestral, reconstruction, sort, convert, full, transitions, recommendation, play, nice, etc, answer, right, worries, couldn, sus, own, mvbm, unfortunately, responsive, inquiries, sorry, yesterday, user, contacted, subject, text, getting, style, wednesday, imo, admittedly, fine, generally, affirms, ought, cited, becomes, tends, affect, measurement, underestimated, overestimated, astonishing, convincing, 9220, 2096, 1930, 4588, 622, 5731, 1219, 2495, 4753, 6045, 5548, 490, 4077, 6367, 8417, 4181, 2777, 6418, 6047, 4362, 4288, 2951, 6848, 2549, 5098, 0866, 4541, 1942, 7713, 2389, 4734, 0244, 4641, 2102, 7708, 1732, 4698, 0304, 4678, 2093, 7854, 1688, 4622, 4686, 2081, 732, 7824, particular, downwards, 0776599, 7045410, 2024869, 7556554, 7203710, 349, 353, 342, 359, 341, summarize, error, estimating, consequence, evident, think, arisen, widely, did, 499556, 500444, 945824, expected, change, gave, accurate, particularly, regard, 632, 1410, 4069, 5262, 0369, 378, 6136, 3327, 7524, 3626, 0282, 258, 4337, 0864, 8329, 5574, 2065, 0202, 247, 0498, 1137, 9956, 5351, 2351, 0203, 245, 0731, 0581, 8858, 6185, 2743, 0387, 242, 3865, 0999, 9000, 6364, 2739, 0330, 9284, 1005, 9045, 6323, 2720, 666, 6320, 14104895, 40690078, 52619893, 03687111, predict, 06588211, 91247554, 49325269, 40862951, 208, 88237689, shows, 182, 194, 224, 206, 209, 216, 214, 192, 186, 232, 178, 213, 180, 228, 210, 239, 207, 196, created, easier, pattern, socket, host, localhost, detectcores, speed, notice, immediately, resemble, 32705, 67295, 278772, imagine, unobserved, terrific, assume, where, comes, theta_a, theta_b, conditioning, depends, total, depth, intend, engineer, simple, circumstances, words, uncertain, among, magnitude, typical, huelsenbeck, pointed, heterogeneous, tend, knowledge, studied, explicitly, talking, documented, couple, entry, princeton, pages,
Text of the page (random words):
alpha alpha setnames rep 0 6 k letters 1 k alpha a b 0 6 0 6 set sigma squared sig2 setnames rep 0 3 k letters 1 k sig2 a b 0 3 0 3 now let s set theta_a and theta_b i m going to use values of 0 and 5 for the two different discrete character states set theta theta setnames c 0 5 letters 1 k theta a b 0 5 now i ll simulate continuous trait data under the multi regime model using our true discrete character history with phytools multiou i m going to assume the starting value is at theta for the root regime this is where the a0 theta y0 comes from simulate continuous trait low x multiou low sim_tree alpha sig2 theta a0 theta y0 head low x t16 t21 t123 t124 t38 t49 0 6548297 0 7918578 0 1488759 0 4502061 0 3959087 0 2302126 terrific so far so good now we re about to imagine that our discrete character history is unobserved so let s pull of just the tip values of our discrete trait from low sim_tree as follows low y as factor getstates low sim_tree type tips head low y t16 t21 t123 t124 t38 t49 a a a a a a levels a b next let s generate a set of stochastic character histories for our discrete state using phytools simmap low mk_fit fitmk phy low y model er pi fitzjohn low mk_fit object of class fitmk fitted or set value of q a b a 0 036172 0 036172 b 0 036172 0 036172 fitted or set value of pi a b 0 32705 0 67295 due to treating the root prior as a nuisance log likelihood 76 278772 optimization method used was nlminb r thinks it has found the ml solution low smps simmap low mk_fit low smps 100 phylogenetic trees with mapped discrete characters let s plot a few of these par mfrow c 5 5 nulo sapply low smps 1 25 plot lwd 1 col cols ftype off direction upwards something that readers will probably notice immediately is that all of these discrete character histories fairly closely resemble each other the true history of our trait next let s fit ouwie models to each of them that ll be fun to speed it along i m going to parallelize across stochastic map trees using the foreach package make our ouwie data frame low ouwie_data data frame genus_species names low x reg low y x low x head low ouwie_data genus_species reg x t16 t16 a 0 6548297 t21 t21 a 0 7918578 t123 t123 a 0 1488759 t124 t124 a 0 4502061 t38 t38 a 0 3959087 t49 t49 a 0 2302126 load foreach and doparallel library foreach library doparallel set up our cluster ncores detectcores 2 ncores 1 14 mc makecluster ncores type psock mc socket cluster with 14 nodes on host localhost registerdoparallel cl mc optimize across all 100 stochastic maps low ouwie_fits foreach i 1 length low smps dopar ouwie ouwie low smps i low ouwie_data model oum simmap tree true root station false stopcluster mc the object we ve created is a long list of fitted ouwie models we can pull out just the specific results that we re interested in to make it easier to see the general pattern foo function x setnames c x theta 1 x solution 1 x loglik c the a the b alpha sigsq log l low ouwie_results t sapply low ouwie_fits foo rownames low ouwie_results 1 nrow low ouwie_results options scipen 5 just for printing round low ouwie_results 1 25 digits 2 the a the b alpha sigsq log l 1 0 08 4 84 0 66 0 42 182 89 2 0 11 4 95 0 49 0 36 194 30 3 0 19 5 07 0 41 0 41 224 99 4 0 20 4 95 0 63 0 48 206 92 5 0 00 4 75 0 67 0 52 209 01 6 0 11 4 95 0 44 0 40 216 47 7 0 08 4 87 0 47 0 41 214 46 8 0 05 4 86 0 48 0 36 192 84 9 0 01 4 91 0 40 0 41 227 67 10 0 03 4 82 0 62 0 41 186 94 11 0 07 4 89 0 50 0 36 193 14 12 0 10 4 86 0 33 0 38 232 95 13 0 03 4 86 0 59 0 37 178 45 14 0 03 4 78 0 43 0 38 213 53 15 0 03 4 88 0 48 0 31 175 21 16 0 10 4 97 0 50 0 33 180 42 17 0 17 4 95 0 49 0 48 228 08 18 0 12 4 92 0 60 0 48 210 78 19 0 04 4 82 0 48 0 37 200 46 20 0 07 4 78 0 42 0 34 199 99 21 0 21 5 29 0 34 0 41 239 15 22 0 05 4 87 0 57 0 40 193 95 23 0 14 5 06 0 47 0 39 207 78 24 0 22 5 21 0 37 0 39 227 25 25 0 10 4 90 0 61 0 43 196 33 options scipen 0 this is pretty cool because it shows that all of our fitted models are pretty close to the generating conditions of mathbf theta 0 5 colmeans low ouwie_results the a the b alpha sigsq log l 0 06588211 4 91247554 0 49325269 0 40862951 208 88237689 we can also compare this to phytools fitmultiou which i predict will get a similar model init setnames c mean low x low y levels low y 1 mean low x low y levels low y 2 log 2 max nodeheights phy var low x max nodeheights phy fitmk phy low y model er rates c theta a theta b alpha sigsq q 1 init theta a theta b alpha sigsq q 1 0 14104895 4 40690078 0 06931472 0 52619893 0 03687111 low fit_mou fitmultiou phy low x low y model er levs 100 parallel true ncores ncores root mle trace 1 maxit 2000 init init iter the a the b alpha sigsq q 1 log l 0 0 1410 4 4069 0 0693 0 5262 0 0369 378 6136 100 0 3327 4 7524 0 3626 0 1910 0 0282 258 4337 200 0 0864 4 8329 0 5574 0 2065 0 0202 247 0498 300 0 1137 4 9956 0 5351 0 2351 0 0203 245 0731 400 0 0581 4 8858 0 6185 0 2743 0 0387 242 3865 500 0 0999 4 9000 0 6364 0 2739 0 0330 241 9284 600 0 1005 4 9045 0 6323 0 2720 0 0331 241 9258 666 0 1009 4 9043 0 6320 0 2719 0 0331 241 9258 done optimizing low fit_mou object of class fitmultiou based on a discretization with k 100 levels fitted multi theta ou model parameters levels a b theta 0 1009 4 9043 alpha 0 632 sigsq 0 2719 estimated q matrix a b a 0 03312624 0 03312624 b 0 03312624 0 03312624 log likelihood 241 9258 r thinks it has found the ml solution even though both analyses gave as parameter estimates quite close to the generating values joint estimation seems to be even more accurate particularly with regard to alpha and sigma 2 that s cool now let s consider the case of a high rate of transition for our discrete character q 0 8 k 2 high q matrix q k k dimnames list letters 1 k letters 1 k diag high q 0 diag high q rowsums high q high q a b a 0 8 0 8 b 0 8 0 8 now the discrete character is expected to change a lot in the true history high sim_tree sim history phy high q anc y0 done simulation s high sim_tree phylogenetic tree with 250 tips and 249 internal nodes tip labels t16 t21 t123 t124 t38 t49 the tree includes a mapped 2 state discrete character with states a b rooted includes branch lengths plot this tree just for fun cols setnames hcl colors n k letters 1 k plot high sim_tree cols ftype off lwd 1 direction upwards ylim c 1 10 par lend 1 legend bottomleft letters 1 k lwd 4 col hcl colors n k cex 0 7 bty n let s generate continuous trait data under the same alpha sigma 2 and mathbf theta 0 5 as before but with our new high q history simulate continuous trait high x multiou high sim_tree alpha sig2 theta a0 theta y0 head high x t16 t21 t123 t124 t38 t49 1 142639 1 634035 1 209871 1 537554 2 871243 1 072081 of course now we might expect our stochastic character histories to differ more from each other and from the true history high y as factor getstates high sim_tree type tips head high y t16 t21 t123 t124 t38 t49 a a a a a a levels a b high mk_fit fitmk phy high y model er pi fitzjohn high mk_fit object of class fitmk fitted or set value of q a b a 0 458785 0 458785 b 0 458785 0 458785 fitted or set value of pi a b 0 499556 0 500444 due to treating the root prior as a nuisance log likelihood 159 945824 optimization method used was nlminb r thinks it has found the ml solution high smps simmap high mk_fit high smps 100 phylogenetic trees with mapped discrete characters much as we did earlier let s plot a few of these just to see what we ve got par mfrow c 5 5 nulo sapply high smps 1 25 plot lwd 1 col cols ftype off direction upwards it should be pretty evident i think that even though they seem to have arisen under the same process in this case the specific details of our discrete character histories vary quite widely one from the other as well as from our generating history i expect that the consequence of this will be that each one of our stochastic character histories is likely to add error and perhaps bias to the estimating of the continuous trait multi regime ou process let s see if that s true make our ouwie data frame high ouwie_data data frame genus_species names high x reg high y x high x head high ouwie_data genus_species reg x t16 t16 a 1 142639 t21 t21 a 1 634035 t123 t123 a 1 209871 t124 t124 a 1 537554 t38 t38 a 2 871243 t49 t49 a 1 072081 once again we ll parallelize across maps using foreach foreach mc makecluster ncores type psock registerdoparallel cl mc optimize across all 100 stochastic maps high ouwie_fits foreach i 1 length high smps dopar ouwie ouwie high smps i high ouwie_data model oum simmap tree true root station false stopcluster mc let s summarize our results foo function x setnames c x theta 1 x solution 1 x loglik c the a the b alpha sigsq log l high ouwie_results t sapply high ouwie_fits foo rownames high ouwie_results 1 nrow high ouwie_results here s just the first 25 again options scipen 5 just for printing round high ouwie_results 1 25 digits 2 the a the b alpha sigsq log l 1 0 67 3 66 0 20 0 73 350 92 2 0 64 3 55 0 19 0 74 352 68 3 0 55 3 73 0 19 0 71 349 64 4 1 11 3 40 0 19 0 75 355 78 5 0 59 3 63 0 19 0 72 350 10 6 1 24 3 54 0 22 0 80 354 80 7 1 51 3 03 0 19 0 78 360 86 8 1 02 3 40 0 21 0 78 354 09 9 1 37 4 09 0 21 0 76 350 52 10 0 94 3 56 0 20 0 75 353 15 11 0 94 4 29 0 21 0 71 342 30 12 1 14 3 25 0 22 0 81 356 49 13 1 66 3 62 0 17 0 74 359 53 14 1 46 3 64 0 20 0 77 356 67 15 1 02 3 13 0 16 0 72 358 30 16 0 31 3 84 0 22 0 72 341 27 17 1 54 3 59 0 19 0 77 358 10 18 1 43 4 06 0 18 0 73 355 63 19 1 53 3 73 0 21 0 78 355 44 20 1 35 4 16 0 18 0 72 352 88 21 1 53 3 08 0 18 0 77 360 42 22 0 42 3 49 0 19 0 73 350 55 23 1 42 3 85 0 21 0 79 354 72 24 0 23 3 70 0 17 0 68 350 52 25 1 02 3 36 0 20 0 77 355 54 options scipen 0 now let s get an average across trees colmeans high ouwie_results the a the b alpha sigsq log l 1 0776599 3 7045410 0 2024869 0 7556554 352 7203710 as we predicted the estimates of theta in particular are biased towards each other but we also see that alpha is biased downwards and sigma 2 upwards finally let s compare this to phytools fitmultiou init setnames c mean high x high y levels high y 1 mean high x high y levels high y 2 log 2 max nodeheights phy var high x max nodeheights phy fitmk phy high y model er rates c theta a theta b alpha sigsq q 1 high fit_mou fitmultiou phy high x high y model er levs 100 parallel true ncores ncores root mle trace 1 maxit 2000 init init iter the a the b alpha sigsq q 1 log l 0 1 9220 3 2096 0 0693 0 1930 0 4588 622 5731 100 1 1219 4 2495 0 4753 0 6045 0 5548 490 4077 200 0 6367 4 8417 0 4181 0 2777 0 6418 473 6047 300 0 4362 4 8024 0 4288 0 2951 0 6848 473 2549 400 0 5098 5 0866 0 4541 0 1942 0 7713 472 2389 500 0 4734 5 0244 0 4641 0 2102 0 7708 472 1732 600 0 4698 5 0304 0 4678 0 2093 0 7854 472 1688 700 0 4622 5 0253 0 4686 0 2081 0 7829 472 1673 732 0 4625 5 0253 0 4682 0 2082 0 7824 472 1673 done optimizing high fit_mou object of class fitmultiou based on a discretization with k 100 levels fitted multi theta ou model parameters levels a b theta 0 4625 5 0253 alpha 0 4682 sigsq 0 2082 estimated q matrix a b a 0 7823842 0 7823842 b 0 7823842 0 7823842 log likelihood 472 1673 r thinks it has found the ml solution this is pretty astonishing seems to make quite a convincing case for joint estimation admittedly this is a very small experiment of course however it does seem to confirm the thesis that if our discrete character changes infrequently then the traditional two step process of generating stochastic character maps then fitting a fixed regime model to each map should work just fine on the other hand it also generally affirms revell 2013 which perhaps ought to be a bit better cited that as the rate of transition in our character goes up the two step process becomes biased and tends to underestimate the difference between the regimes this also may affect the measurement of alpha which is underestimated and sigma 2 which is overestimated this bias goes away when we jointly model x and y using our new finite space approximation which is very cool imo ok more on this later posted by liam revell at 3 42 pm no comments email this blogthis share to x share to facebook share to pinterest wednesday july 29 2026 stochastic character mapping under the threshold model using fitthresh in phytools unfortunately i m not usually so responsive to email inquiries these days sorry but just yesterday a phytools user contacted me with the subject line getting a simmap style tree from a threshold model result and the following email text i have a weird question regarding your threshold model functions i have a continuous trait i am making discrete bins out of to test whether multivariate rates of shape evolution vary with the levels of this now discrete trait using mvmorph fitthresh makes a lot of sense for reconstructing this character and results in ancestral values that look different more sensible than make simmap based on a reconstruction of the underlying continuous character but i need the threshold model results in some sort of format that will work with mvbm e g simmap is there a not too complicated way to convert the full q matrix into transitions between the original discrete character or do you have a better recommendation for making fitthresh play nice with the trait fitting functions in phytools mvmorph ouwie etc if there is no good answer right now no worries just couldn t sus it out on my own indeed this is something that i d thought of already that is that the discrete diffusion finite space approximation that we use in fitthresh and other new methods of phytools would enable new sorts of stochastic character mapping including of the threshold model and indeed of actual continuous traits before i go on to show how to do this i will note that this is an approximation of the corresponding stochastic character map for the threshold process only inasmuch as technically the threshold is crossed infty times by the underlying fractal bm process when it moves from one side of the threshold to the other if you didn t already know this just trust me it s true so how do we do this well so far i have not yet automated this in software but it s pretty straightforward if we go through it step by step i m going to start with the simplest example which is the threshold model with two states let s call them a and b load phytools library phytools first let s simulate some data under this process simulate tree phy pbtree n 80 scale 1 phy phylogenetic tree with 80 tips and 79 internal nodes tip labels t28 t31 t32 t23 t24 t19 rooted includes branch length s simulate liabilities liability fastbm phy a 0 5 head liability 10 t28 t31 t32 t23 t24 t19 t25 t26 0 9810162 1 1667724 0 9394853 0 5019343 0 2609373 0 1475927 0 4670823 0 2442964 t27 t33 0 7697432 1 0546414 threshold liabilities thresh function x if x 0 a else b x sapply liability thresh head x 20 t28 t31 t32 t23 t24 t19 t25 t26 t27 t33 t35 t37 t38 t79 t80 t9 t29 t30 t46 t47 b b b b b a a b a a a a a a a b a a a a now let s proceed to fit ...
|
