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):
r dependent continuous trait evolution under a multi regime ou model a couple of days ago on this blog i posted an entry on comparing a discrete character dependent multi regime joint ou model to ouwie when the regime history is known for those not following closely at home the joint discrete continuous trait model that i m talking about uses a finite space discretized diffusion approximation as documented in prior posts to this blog e g 1 2 3 4 5 following boucher démery 2016 also see our in review biorxiv pre print a more typical workflow with ouwie and other similar regime based methods is to generate a set of stochastic character maps of the discrete character regime following huelsenbeck et al 2003 fit a multi regime ou model to each tree in the set and then average the results across maps back in 2013 i pointed out that for rate heterogeneous models this would tend to result in an underestimate of the difference in rates sigma 2 between regimes the same seems very likely to be true for multi regime ou models but to my knowledge this has not been studied explicitly we might hope or expect that this bias would go away if we can jointly model the discrete continuous characters the magnitude of this bias is also likely to vary as a function of the transition rate for the discrete trait in other words the higher the transition rates in q of our discrete trait the more uncertain our discrete character history and the more it will vary among stochastic character maps and the more each map will differ from the unknown true history of our trait on a sunday morning i thought it would be interesting to engineer a simple comparison between this traditional ouwie workflow and fitmultiou both for circumstances in which the transition rate of the discrete character is low high to begin with let s load the packages we intend to use load packages library phytools library ouwie now i m going to simulate a tree a discrete character history on that tree with a low rate keeping in mind of course that there is no special numeric value of q that is low or high it all depends on the total depth of our tree phy pbtree n 250 scale 10 phy phylogenetic tree with 250 tips and 249 internal nodes tip labels t16 t21 t123 t124 t38 t49 rooted includes branch length s q 0 05 k 2 low q matrix q k k dimnames list letters 1 k letters 1 k diag low q 0 diag low q rowsums low q low q a b a 0 05 0 05 b 0 05 0 05 low sim_tree sim history phy low q anc y0 sample letters 1 k 1 done simulation s low 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 low 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 great clearly we ve simulated a very low rate of discrete character evolution on this tree now let s simulate a multi regime ou process conditioning on this true discrete character history to do this we ll use multiou in phytools even though we re using the same value of alpha and sigma 2 for both of our discrete character regimes we need to set values of these two parameters for each of the two character levels set 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...
|
