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):
rmost threshold is fixed totally just for fun let s fit compare a symmetric m k model to see if they re distinguishable based on these data mk_model fitmk phy x model er anova mk_model thresh_fit log l d f aic weight mk_model 36 04421 1 74 08842 0 1050167 thresh_fit 33 90152 1 69 80305 0 8949833 now to undertake stochastic mapping using our threshold model we need to pull out the implicit m k model that s hidden within our fitthresh object mkm thresh_fit mk_fit i m not going to print it because the q matrix is 200 times 200 but readers following along should feel free to go for it next because this hidden object is not quite the same as a standard fitmk object let s add the element root prior as follows mkm root prior fitzjohn now believe it or not we re totally ready for stochastic mapping when we run the following code we re probably going to see the message warning in rstate p sum p some probabilities slightly 0 setting p 0 to zero this is just because our matrix is so big some of the elements are close to zero so we can safely ignore it i m only going to do one stochastic character map here typically we should do many of these of course this would be accomplished by modifying nsim and then iterating all subsequent steps over each stochastic map tree smp simmap mkm nsim 1 pi mle warning in rstate p sum p some probabilities slightly 0 setting p 0 to zero warning in rstate p sum p some probabilities slightly 0 setting p 0 to zero warning in rstate p sum p some probabilities slightly 0 setting p 0 to zero warning in rstate p sum p some probabilities slightly 0 setting p 0 to zero once again this is a big object because it contains all 200 of the finely discretized levels of our diffusion approximation effectively a stochastic map of liabilities to turn this into a stochastic character map of our original discrete trait we need to simply merge liabilities on each side of the single threshold as follows ss colnames mkm data merged smp mergemappedstates smp ss which as numeric ss 0 a merged smp mergemappedstates merged smp ss which as numeric ss 0 b merged smp phylogenetic tree with 80 tips and 79 internal nodes tip labels t28 t31 t32 t23 t24 t19 the tree includes a mapped 2 state discrete character with states a b rooted includes branch lengths great let s plot it plot merged smp direction upwards ftype off lwd 3 colors setnames hcl colors n 2 c a b legend bottomleft c a b col hcl colors n 2 lwd 2 bty n we can see that this is clearly capturing a form of the threshold model history of our threshold character by seeing just how different it looks compared to running the same analysis with a standard m k model mk smp simmap mk_model nsim 1 plot mk smp direction upwards ftype off lwd 3 colors setnames hcl colors n 2 c a b legend bottomleft c a b col hcl colors n 2 lwd 2 bty n we see for instance that threshold crossing usually involves many switches back forth under the threshold model and none in the m k model which makes perfect sense things get only a little more complicated when i want to fit the multi state model let s see if i can run through that first i ll reuse my original tree liabilities to simulate some tip data as follows y threshstate liability setnames c 0 1 2 letters 1 3 head y t28 t31 t32 t23 t24 t19 b c b b b a next fit the multi state threshold model as follows ms_thresh fitthresh phy y sequence letters 1 3 ms_thresh object of class fitthresh set value of sigsq of the liability 1 0 set or estimated threshold s 0 7568 0 341537 log likelihood 51 847661 lowermost threshold is fixed note that the lower estimated threshold is not a separately estimable parameter otherwise the model would be non identifiable it s normally set to zero but fitthresh does something different that i ll probably fix in the future it centers the whole liability distribution on zero instead since the liabilities are unitless and scaleless this doesn t matter at all for the model fit or relative positions of the thresholds but our results are a bit more annoying to interpret pull out our hidden m k object do stochastic mapping as before mkn ms_thresh mk_fit mkn root prior fitzjohn smp simmap mkn nsim 1 warning in rstate p sum p some probabilities slightly 0 setting p 0 to zero warning in rstate p sum p some probabilities slightly 0 setting p 0 to zero warning in rstate p sum p some probabilities slightly 0 setting p 0 to zero warning in rstate p sum p some probabilities slightly 0 setting p 0 to zero now repeat the mergemappedstates step but this time using the estimated thresholds ss mkn states merged smp 2 mergemappedstates smp ss which as numeric ss ms_thresh threshold 1 a merged smp 2 mergemappedstates merged smp 2 ss which as numeric ss ms_thresh threshold 1 as numeric ss ms_thresh threshold 2 b merged smp 2 mergemappedstates merged smp 2 ss which as numeric ss ms_thresh threshold 2 c merged smp 2 phylogenetic tree with 80 tips and 79 internal nodes tip labels t28 t31 t32 t23 t24 t19 the tree includes a mapped 3 state discrete character with states a b c rooted includes branch lengths let s plot it plot merged smp 2 direction upwards ftype off lwd 3 colors setnames hcl colors n 3 c a b c legend bottomleft c a b c col hcl colors n 3 lwd 3 bty n wow that s beautiful very cool posted by liam revell at 2 23 pm no comments email this blogthis share to x share to facebook share to pinterest comparing a discrete character dependent multi regime ou joint model to ouwie in some recent posts to this blog e g 1 2 3 4 i have described a new discrete character dependent multi optimum model in phytools that uses the finite space or discrete diffusion approximation of boucher démery 2016 also see our in review biorxiv pre print even though they implement different models i thought it might be interesting to compare this new method fitmultiou to a fixed regime multi optimum ou model fit using the popular ouwie package by jeremy beaulieu and brian o meara since our new method integrates over uncertainty in the regime history by jointly optimizing an m k transition process for the regimes along with the multi optimum stochastic process of our continuous trait my thesis is that the parameter estimates we obtain should be highly similar between the two models if our discrete trait changes infrequently to explore this we can start by simulating a tree a discrete character history and some data using phytools as follows load phytools library phytools simulate a tree n 100 number of taxa phy pbtree n n scale 10 phy phylogenetic tree with 100 tips and 99 internal nodes tip labels t2 t3 t89 t90 t35 t53 rooted includes branch length s set the transition matrix of our discrete trait q 0 05 k 2 q matrix q k k dimnames list letters 1 k letters 1 k diag q 0 diag q rowsums q q a b a 0 05 0 05 b 0 05 0 05 simulate trait history sim_tree sim history phy q anc sample letters 1 k 1 done simulation s cols setnames hcl colors n k letters 1 k plot sim_tree cols ftype off lwd 2 direction upwards par lend 1 legend bottomleft letters 1 k lwd 4 col hcl colors n k cex 0 7 bty n this is great because our tree has very few changes in the character my hypothesis is that this will make parameter estimates of the ou process very similar between a fixed regime model as in ouwie and our new fitmultiou method now let s set the generating conditions of our discrete character dependentn multi theta ornstein uhlenbeck process the parameters alpha and sigma 2 will be the same across all regimes though we still have to specify k 2 of each for our generator while theta will vary according to the state of our discrete trait set alpha alpha setnames rep 0 3 k letters 1 k alpha a b 0 3 0 3 set sigma squared sig2 setnames rep 0 1 k letters 1 k sig2 a b 0 1 0 1 set theta theta setnames c 0 5 2 letters 1 k theta a b 0 5 2 0 at this point we can simulate our continuous trait using phytools multiou simulate continuous trait x multiou sim_tree alpha sig2 theta a0 0 head x t2 t3 t89 t90 t35 t53 0 3416784 0 6347880 0 3128115 0 5001694 0 7747887 0 5183583 our discrete character is already simulated but we need to pull of the tip values from the sim_tree simmap object to use them in our analysis pull off discrete trait y as factor getstates sim_tree tips head y t2 t3 t89 t90 t35 t53 a a a a a a levels a b for completeness let s start by fitting our null model using fitmultiou we can then confirm that this fitted model matches to a reasonable degree they will only converge exactly as levs goes towards infty what we d obtain using geiger fitcontinuous and phytools fitmk under the same model assumptions i m going to set levs 200 for this analysis but this actually takes a very long time to run much more than twice as long as levs 100 fit null model fit_null fitmultiou phy x y model er levs 200 parallel true ncores 10 root mle trace 1 null_model true iter theta alpha sigsq q 1 log l 0 1 9456 0 0721 1 1329 0 0914 175 6004 100 0 6244 0 0116 0 0897 0 0271 98 3007 200 0 6123 0 0053 0 0892 0 0283 98 2222 271 0 6338 0 0004 0 0860 0 0292 98 1838 done optimizing fit_null object of class fitmultiou based on a discretization with k 200 levels fitted multi theta ou model parameters levels a b theta 0 6338 alpha 4e 04 sigsq 0 086 estimated q matrix a b a 0 02921177 0 02921177 b 0 02921177 0 02921177 log likelihood 98 1838 r thinks optimization may not have converged fit single regime ou model using fitcontinuous ou_fit geiger fitcontinuous phy x model ou ou_fit now again let s compare to geiger fitcontinuous and phytools fitmk geiger fitted comparative model of continuous data fitted ou model parameters alpha 0 000000 sigsq 0 085106 z0 0 625022 model summary log likelihood 68 250476 aic 142 500953 aicc 142 750953 free parameters 3 convergence diagnostics optimization iterations 100 failed iterations 0 number of iterations with same best fit 50 frequency of best fit 0 500 object summary lik likelihood function bnd bounds for likelihood search res optimization iteration summary opt maximum likelihood parameter estimates fit mk model using fitmk mk_fit fitmk phy y model er pi equal mk_fit object of class fitmk fitted or set value of q a b a 0 028064 0 028064 b 0 028064 0 028064 fitted or set value of pi a b 0 5 0 5 due to treating the root prior as a flat log likelihood 29 168076 optimization method used was nlminb r thinks it has found the ml solution you can compare the model parameter estimates but let s also assure ourselves that the likelihood seems to be converging on the same value as follows compute null log l from fitcontinuous fitmk results null_logl loglik ou_fit loglik mk_fit attr null_logl df 4 fix d f null_logl 1 97 41855 attr df 1 4 compare to fitmultiou loglik fit_null 1 98 18378 attr df 1 4 this is pretty close again we would expect these two values to get even closer for higher levs but as currently implemented this already takes a really long time to run ok now let s bring ouwie into the picture we can start by loading the package which i recently updated from cran load ouwie library ouwie now for ouwie we need to put our data in a special data frame format as follows compile our data for ouwie ouwie data data frame genus_species names x reg y x x head ouwie data genus_species reg x t2 t2 a 0 3416784 t3 t3 a 0 6347880 t89 t89 a 0 3128115 t90 t90 a 0 5001694 t35 t35 a 0 7747887 t53 t53 a 0 5183583 why don t we start by simply re fitting our null ou model in ouwie this should give us a result that quite closely matches what we obtained using geiger fitcontinuous i ll still give it our known discrete character data history but i set model ou1 to specify that i want a single theta model only fit ouwie null model fitou smp ouwie sim_tree ouwie data model ou1 simmap tree true root station false warning an algorithm was not specified defaulting to computing the determinant and inversion of the vcv initializing finished begin thorough search finished summarizing results fitou smp fit lnl aic aicc bic model ntax 68 25048 142 501 142 751 150 3165 ou1 100 rates alpha sigma sq 1 524044e 08 8 510672e 02 optima 1 estimate 0 6250224 se 0 3348563 half life another way of reporting alpha alpha 45480781 arrived at a reliable solution hopefully we see that this fitted model pretty closely matches what we got using fitcontinuous earlier next i m going to go ahead fit our discrete character dependent multi theta ou model using phytools fitmultiou this is the model that i ve been blogging about recently but to remind the reader this is a joint discrete continuous trait model not the fixed regime model of ouwie but i m hypothesizing that our continuous trait model parameter estimates should pretty closely match what we d get from ouwie using the true history or for that matter a stochastic character history just because our discrete character changes so infrequently on the tree to start with i m going to try to get reasonable starting values for my fitmultiou parameters as follows identify sensible starting parameter values init setnames c mean x y levels y 1 mean x y levels y 2 log 2 max nodeheights phy var x max nodeheights phy fitmk phy y model er rates c theta a theta b alpha sigsq q 1 init theta a theta b alpha sigsq q 1 0 27841916 1 61842864 0 06931472 0 12828760 0 02806402 then i can go ahead and fit the joint model fit discrete trait dependent model fit_mou fitmultiou phy x y model er levs 200 parallel true ncores 10 root mle trace 1 maxit 2000 init init iter the a the b alpha sigsq q 1 log l 0 0 2784 1 6184 0 0693 0 1283 0 0281 92 1417 100 0 4558 1 8575 0 3260 0 1157 0 0268 72 4242 200 0 4520 1 8537 0 3254 0 1143 0 0267 72 4172 300 0 4531 1 8535 0 3252 0 1145 0 0269 72 4170 400 0 4520 1 8523 0 3246 0 1142 0 0271 72 4163 404 0 4519 1 8523 0 3246 0 1141 0 0271 72 4163 done optimizing fit_mou object of class fitmultiou based on a discretization with k 200 levels fitted multi theta ou model parameters levels a b theta 0 4519 1 8523 alpha 0 3246 sigsq 0 1141 estimated q matrix a b a 0 02713072 0 02713072 b 0 02713072 0 02713072 log likelihood 72 4163 r thinks it has found the ml solution finally we can fit our fixed regime model using the ouwie package here don t forget that though i will be using the true discrete character history this is almost never known in practice indeed if we genuinely knew the discrete character history i would always recommend fitting a fixed regime model not a joint model fit multi regime ou model using ouwie fitoum smp ouwie sim_tree ouwie data model oum simmap tree true root station false warning an algorithm was not specified defaulting to computing the determinant and inversion of the vcv initializing finished begin thorough search finished summarizing results fitoum smp fit lnl aic aicc bic model ntax 44 47594 96 95187 97 37292 107 3726 oum 100 rates b a alpha 0 3007500 0 3007500 sigma sq 0 1198701 0 1198701 optima b a estimate 1 75307314 0 5184579 se 0 09854245 0 1132707 half life another way of reporting alpha b a 2 304729 2 304729 arrived at a reliable solution cool now if w...
|
