[med-svn] [r-cran-ape] 01/02: Imported Upstream version 3.3
Dylan Aïssi
bob.dybian-guest at moszumanska.debian.org
Thu Jun 18 17:15:40 UTC 2015
This is an automated email from the git hooks/post-receive script.
bob.dybian-guest pushed a commit to branch master
in repository r-cran-ape.
commit c682bf577497b2df37d9580c32c044796325bb2d
Author: Dylan Aïssi <bob.dybian at gmail.com>
Date: Thu Jun 18 19:14:37 2015 +0200
Imported Upstream version 3.3
---
DESCRIPTION | 14 +-
MD5 | 98 ++++++------
NAMESPACE | 5 +
R/CDF.birth.death.R | 174 +++++++++++++++++++++-
R/DNA.R | 29 ++--
R/PGLS.R | 7 +-
R/ace.R | 14 +-
R/as.bitsplits.R | 10 +-
R/as.phylo.R | 7 +-
R/binaryPGLMM.R | 293 ++++++++++++++++++++++++++++++++++++
R/compar.gee.R | 7 +-
R/corphylo.R | 355 ++++++++++++++++++++++++++++++++++++++++++++
R/dbd.R | 9 +-
R/dist.topo.R | 72 +++++----
R/parafit.R | 6 +-
R/plot.phylo.R | 38 +++--
R/read.tree.R | 8 +-
R/root.R | 8 +-
R/write.nexus.R | 15 +-
R/write.nexus.data.R | 6 +-
build/vignette.rds | Bin 194 -> 192 bytes
data/hivtree.newick.rda | Bin 2190 -> 2189 bytes
data/landplants.newick.rda | Bin 576 -> 574 bytes
data/opsin.newick.rda | Bin 446 -> 444 bytes
inst/doc/MoranI.pdf | Bin 237986 -> 238306 bytes
man/DNAbin.Rd | 2 +-
man/LTT.Rd | 6 +-
man/MoranI.Rd | 4 +-
man/ace.Rd | 14 +-
man/as.bitsplits.Rd | 8 +-
man/binaryPGLMM.Rd | 268 +++++++++++++++++++++++++++++++++
man/boot.phylo.Rd | 3 +-
man/coalescent.intervals.Rd | 7 +-
man/collapsed.intervals.Rd | 13 +-
man/corphylo.Rd | 224 ++++++++++++++++++++++++++++
man/def.Rd | 2 +-
man/del.gaps.Rd | 32 ++--
man/dist.topo.Rd | 18 ++-
man/is.binary.tree.Rd | 2 +-
man/landplants.Rd | 3 +-
man/node.depth.Rd | 10 +-
man/read.dna.Rd | 6 +-
man/read.tree.Rd | 2 +-
man/rlineage.Rd | 53 +++++--
man/root.Rd | 11 ++
man/seg.sites.Rd | 5 -
man/skyline.Rd | 42 ++----
man/skylineplot.Rd | 34 +----
man/where.Rd | 3 +-
man/woodmouse.Rd | 3 +-
man/write.dna.Rd | 6 +-
src/dist_dna.c | 24 ++-
52 files changed, 1653 insertions(+), 327 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index e416714..1549b07 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
Package: ape
-Version: 3.2
-Date: 2014-12-05
+Version: 3.3
+Date: 2015-05-28
Title: Analyses of Phylogenetics and Evolution
Authors at R: c(person("Emmanuel", "Paradis", role = c("aut", "cre", "cph"), email = "Emmanuel.Paradis at ird.fr"),
person("Simon", "Blomberg", role = c("aut", "cph")),
@@ -13,6 +13,7 @@ Authors at R: c(person("Emmanuel", "Paradis", role = c("aut", "cre", "cph"), email
person("Julien", "Dutheil", role = c("aut", "cph")),
person("Olivier", "Gascuel", role = c("aut", "cph")),
person("Christoph", "Heibl", role = c("aut", "cph")),
+ person("Anthony", "Ives", role = c("aut", "cph")),
person("Daniel", "Lawson", role = c("aut", "cph")),
person("Vincent", "Lefort", role = c("aut", "cph")),
person("Pierre", "Legendre", role = c("aut", "cph")),
@@ -29,10 +30,11 @@ Depends: R (>= 3.0.0)
Suggests: gee, expm
Imports: nlme, lattice, graphics, stats, tools, utils, parallel
ZipData: no
-Description: ape provides functions for reading, writing, plotting, and manipulating phylogenetic trees, analyses of comparative data in a phylogenetic framework, ancestral character analyses, analyses of diversification and macroevolution, computing distances from allelic and nucleotide data, reading and writing nucleotide sequences as well as importing from BioConductor, and several tools such as Mantel's test, generalized skyline plots, graphical exploration of phylogenetic data (alex [...]
+Description: Functions for reading, writing, plotting, and manipulating phylogenetic trees, analyses of comparative data in a phylogenetic framework, ancestral character analyses, analyses of diversification and macroevolution, computing distances from allelic and nucleotide data, reading and writing nucleotide sequences as well as importing from BioConductor, and several tools such as Mantel's test, generalized skyline plots, graphical exploration of phylogenetic data (alex, trex, krono [...]
License: GPL (>= 2)
URL: http://ape-package.ird.fr/
-Packaged: 2014-12-04 21:10:22 UTC; paradis
+NeedsCompilation: yes
+Packaged: 2015-05-29 09:20:40 UTC; paradis
Author: Emmanuel Paradis [aut, cre, cph],
Simon Blomberg [aut, cph],
Ben Bolker [aut, cph],
@@ -44,6 +46,7 @@ Author: Emmanuel Paradis [aut, cre, cph],
Julien Dutheil [aut, cph],
Olivier Gascuel [aut, cph],
Christoph Heibl [aut, cph],
+ Anthony Ives [aut, cph],
Daniel Lawson [aut, cph],
Vincent Lefort [aut, cph],
Pierre Legendre [aut, cph],
@@ -57,6 +60,5 @@ Author: Emmanuel Paradis [aut, cre, cph],
Korbinian Strimmer [aut, cph],
Damien de Vienne [aut, cph]
Maintainer: Emmanuel Paradis <Emmanuel.Paradis at ird.fr>
-NeedsCompilation: yes
Repository: CRAN
-Date/Publication: 2014-12-05 17:18:03
+Date/Publication: 2015-05-29 12:44:18
diff --git a/MD5 b/MD5
index 70c2063..d6d5cdd 100644
--- a/MD5
+++ b/MD5
@@ -1,26 +1,27 @@
eb723b61539feef013de476e68b5c50a *COPYING
-85f4448ef7bc83a44650c96e66d1eec8 *DESCRIPTION
-c529a081c74a3eb99da54b438ee7f2b4 *NAMESPACE
+50b68c64f7094784d90ab4c5a2e7d44a *DESCRIPTION
+ce33ce7a078a76ba0b7d15c174e19078 *NAMESPACE
854a025cb7e5da3e4fe230c0be950d08 *NEWS
0c7bc9101516fd26fb3ddbedbe25b6a3 *R/CADM.global.R
e042efca1d8a00d4178204f5f481b64d *R/CADM.post.R
-5fd68920346c5aca7e7695fd9dbc1c3c *R/CDF.birth.death.R
+f6dfdae2ca8529a2f0ba1de6b90fed0f *R/CDF.birth.death.R
fdb9cfa0cbda82bda982b290693e44e3 *R/Cheverud.R
-94ab10b6cd1040fb0c4937da0f5f1a50 *R/DNA.R
+144ad46a2f5fcc069f97971f4d405b91 *R/DNA.R
0fbc7715dfdc8d54ded16d78ca3100f8 *R/MPR.R
bb95af56d882b6162aa517a77140175e *R/MoranI.R
-dd2d159ad60c382ecc7331ddb10b3e5b *R/PGLS.R
+218b8cf3c13a700c757b7b2303dc897e *R/PGLS.R
a683ed4e25ad199edc7d30e7b55181a6 *R/SDM.R
6b54088bf2514e3aa70edf05fe88127d *R/SlowinskiGuyer.R
-a06252b90c1e32e72619858225f0b51d *R/ace.R
+eaad016718d9ee8c79b15d252879072c *R/ace.R
4ce79cf3f3ff49bef989454d86d0c891 *R/additive.R
9fca35b1005cce7b2a260b553dd971a3 *R/alex.R
9fe874382f024a98f62a0ccfcd6d09ac *R/all.equal.phylo.R
-73f08934c47e0dfb29c49a9d0885cdb3 *R/as.bitsplits.R
+06b049b681ab2e594124a2b010c5cddd *R/as.bitsplits.R
c94018d5e792c72e20ce84085b2df9e7 *R/as.matching.R
-7989b47a70f6003ac477370d92e9655e *R/as.phylo.R
+91674814efe984a1dcc90b3d7d8f425f *R/as.phylo.R
1bb4dcfb44a7485710ec7f599d4b2ee8 *R/as.phylo.formula.R
29efb767d4a207910c1da99c1eb6cb69 *R/balance.R
+607d7d6bc0ec5559ce1de5b1d1aab0c1 *R/binaryPGLMM.R
76b5ca45b7210c85fdc20cd1579b33a2 *R/bind.tree.R
a28930cdae07a271403e988e86a0c324 *R/biplot.pcoa.R
92d9db6d7c1150d1a7c9443fd3d5cb04 *R/birthdeath.R
@@ -33,17 +34,18 @@ f6f25bb8f94f424dd6bbbab64eedb00d *R/clustal.R
4f52fb37d822be586f0c5971295167dc *R/coalescent.intervals.R
054bec70587f0056636b2cc52754c7f9 *R/collapse.singles.R
338accc0accb8e67f31db5a90c033d2f *R/collapsed.intervals.R
-af76fbb797185c9738442161a2c923b4 *R/compar.gee.R
+01e979242ba4667e48a80fd32c03f254 *R/compar.gee.R
89ce53eb1bb8c0fb36e95810cf7cd769 *R/compar.lynch.R
207154a3a9b9ebbe5d7c29995565dc82 *R/compar.ou.R
8d7c71929156744fd0e4fac370bf9456 *R/compute.brtime.R
0092074917dc5631dc62fb9a3016145c *R/cophenetic.phylo.R
13f2dac84d7b8a1da3df0e0c11b4ab1c *R/cophyloplot.R
-d5ae510a60a962ca6420366e3bc9cb1a *R/dbd.R
+fd39268020a494980293c92d0971caca *R/corphylo.R
+1c3460b48ed1e772c40e5993cd663d43 *R/dbd.R
3822f0bb0a9ed4c8c19654e86ef34359 *R/def.R
93480a5b64e0d37f5450891629557615 *R/delta.plot.R
dfd5bb35f1cb1fd9154d023e0e4cfc2b *R/dist.gene.R
-3d29294945ed629de727745d614cc0dd *R/dist.topo.R
+eba6cb1e8dd0359568a46907c445351e *R/dist.topo.R
b28ced504fedeb7f991f7eba10ad06df *R/diversi.gof.R
8b2ec4004022afdc7e2cb42f2657b628 *R/diversi.time.R
fd662cd4bbd6ce536ef379eca2e516ab *R/drop.tip.R
@@ -74,12 +76,12 @@ eb5dd3843951b8bde1b8d53ab2d1353a *R/multi2di.R
e3f22d0f260c43be87a17b0ab091e2bb *R/njs.R
080c82acd4a90a88a1bbd098ebfe6b60 *R/nodelabels.R
ae2aeb0e8aef7f8d4b19939ca61b3482 *R/nodepath.R
-7425d422f7ca9b3966bb65a4e4336bbe *R/parafit.R
+d624cca14aa569e94e1f1e13e9091337 *R/parafit.R
fc0260f9e8e17fc1446117580bbc93cc *R/pcoa.R
4999c85e82105ded85a214204321c3d1 *R/phydataplot.R
e71db002c66c277bfb57f6914ca143d4 *R/phymltest.R
61d1360195bf4d5e6e2d8c30de335dc8 *R/pic.R
-0b0507e9993d4a68aea85aa4ff1ac07b *R/plot.phylo.R
+9ff3970a2edd690be453c7f06076341b *R/plot.phylo.R
e579ec65c808c29e1ecaae1501784b37 *R/plot.popsize.R
736506e2f67e90cf9326abead437d298 *R/plotPhyloCoor.R
1e2485437566ca9af99d93b4580cbbc2 *R/print.lmorigin.R
@@ -90,10 +92,10 @@ b13dfb8f455b1c9e74a364085f72dbce *R/read.caic.R
c1fc8ab4715f1d58145f533fbaf776e4 *R/read.dna.R
56361b2157d0afb1050233a6ab8ee1b5 *R/read.nexus.R
13ce7f5c7d1bcb7101469d12651e99c8 *R/read.nexus.data.R
-23541fbea792b010f52b89735a660582 *R/read.tree.R
+88e633c80886f5a6c9f2b34918452d15 *R/read.tree.R
86361682a3d477e8e63fff53dc67279e *R/reconstruct.R
b1218f7862a66f566e4f3165f28fe47c *R/reorder.phylo.R
-a71f335ebf44985121b28bfcf86d698d *R/root.R
+d1b097835752cd0c7672a68789bdb610 *R/root.R
f584366b32e7414c669714ba5b84951b *R/rotate.R
b7158b84c7ee7d9dcb2e0abeb6005cb0 *R/rtree.R
0730a0bed794d8297576e732868d23ec *R/rtt.R
@@ -112,14 +114,14 @@ a40ae9ad30c221d4ed14b90e8b406f93 *R/vcv.phylo.R
31b3bb1feed474692f07fcebe3a61ac7 *R/vcv2phylo.R
16d5511c871e41b9b2673052e03c8b33 *R/which.edge.R
aa484437bf64a39b737665b0769e320b *R/write.dna.R
-846c603fdb8bda22e5d3e2c81681221f *R/write.nexus.R
-c2006919805d0628cacc83db2bb7eff0 *R/write.nexus.data.R
+4e5d7a7904fe6496a44c4b69f674290b *R/write.nexus.R
+89496c0df70205654456d3b7bb7ba41f *R/write.nexus.data.R
17d72a136a8131ea46ef9a20fcfd4d36 *R/write.tree.R
a918c086a449bcca5ccbb04f7e6d50a9 *R/yule.R
c8d3aa3fe64e75e61af07a1b11c74f3f *R/yule.time.R
1eb44ff9e5a036eb845faa1598ce5009 *R/zoom.R
3387c0d0c1f913f8471e1bb34bd2e516 *R/zzz.R
-2b7357b75c0beefb393056c925cc168e *build/vignette.rds
+85eb7d45a064be571370be6e598c3cf5 *build/vignette.rds
db9083e8750aff839d5ebf3ed982f1f1 *data/HP.links.rda
9d9f9232839665422709ded1e541d038 *data/bird.families.rda
a14a6df0f3a735ebc056065077788c90 *data/bird.orders.rda
@@ -127,29 +129,29 @@ f74f9ed80c04756021cc093d40ca9ff9 *data/carnivora.csv.gz
4eaf8cbaefa2e8f8d395a9b482ee9967 *data/chiroptera.rda
1c74c3b99d08b0e17eea3ec1065c12d2 *data/cynipids.rda
7fe760c2f3b4deba0554aae6138cb602 *data/gopher.D.rda
-f0fa8d1eab01b77eec3f4ec6552369f7 *data/hivtree.newick.rda
+41a3ee7725bfad5558773b3ef4ad7eb0 *data/hivtree.newick.rda
8d14f95319d0a5cdc8faa60a1d0085ce *data/hivtree.table.txt.gz
-cc483099ae0a0386f4b8c0e4905f15c8 *data/landplants.newick.rda
+20d7287583c3ba5ab039a770d3172fc1 *data/landplants.newick.rda
31be81fe3faca11f98d3e74c090bc59e *data/lice.D.rda
38edbd84a0a067322c40db8d71fb1289 *data/lmorigin.ex1.rda
e3ce9e3444182fea2e65df2e150ea0db *data/lmorigin.ex2.rda
ce7a56faebdf286fdf5ba6c8c3699a79 *data/mat3.RData
e2d1339025ed901009bfed58dc6505ff *data/mat5M3ID.RData
101d0ab2e981b0987cde704a2dee1d8d *data/mat5Mrand.RData
-b0bed24573a2dc201f9e82ab3dba6599 *data/opsin.newick.rda
+34b42702379db0d59db817f7c6c33709 *data/opsin.newick.rda
39e4fece2bdc527d7a9d4d83d023a947 *data/woodmouse.rda
828290996b613493f96e0fab024beebb *inst/CITATION
3f54f3775bcf382e25df2a12228894f6 *inst/doc/MoranI.R
abf2e3e84707634e32f3e7cc6a5ee7ae *inst/doc/MoranI.Rnw
-b2ab7d195a318f9f10f2e3801c73ebfd *inst/doc/MoranI.pdf
+d00d528ac1f7e054d372178a314558de *inst/doc/MoranI.pdf
e6876b193a0df06697c788a8e48cf4bc *man/CADM.global.Rd
-060222d269ef45f0c7df76bb723e7d45 *man/DNAbin.Rd
+b79f29b6efb97d7d4cb5408810a23893 *man/DNAbin.Rd
8b9bc214e32cde4c6e5970e48ff30c5f *man/Initialize.corPhyl.Rd
-4af2b7b7c40ed3ade258a0f08ddb16c0 *man/LTT.Rd
+4f9582ee3a35e609a01779e1b1bf32e7 *man/LTT.Rd
ff05fd6fa0a2750b53abf025cdc021d3 *man/MPR.Rd
-e63ecbd9a4cf6de364e595d5999e3a2e *man/MoranI.Rd
+303135aa8664be5cb518e0cbf2018b2c *man/MoranI.Rd
17486c0fd29fb6f4a752c53fe37142c4 *man/SDM.Rd
-5027a986894645f46bade673bf835a3e *man/ace.Rd
+b4b2d3c708bf1ebca6f1f3b4cb752243 *man/ace.Rd
60adac0f324c7b08a756d8f45c02b024 *man/add.scale.bar.Rd
32d8ad6eeea971901aee33f1f654ae7c *man/additive.Rd
0f37535a86413945bf3205b7a7b7c20b *man/alex.Rd
@@ -157,7 +159,7 @@ e63ecbd9a4cf6de364e595d5999e3a2e *man/MoranI.Rd
416c995c868f02d17fe9dc14a106b71c *man/ape-internal.Rd
587cbd9d9a326b3990ed20bfe5016165 *man/ape-package.Rd
5bba4ae4bfc66b613855cfc182d9b1bc *man/as.alignment.Rd
-73c00aae390eaad1a86f3ede4a41e1a4 *man/as.bitsplits.Rd
+2771f272ad46b279022892d9f5d21eb2 *man/as.bitsplits.Rd
4f014cf2923e2eab6188acd48e8096fa *man/as.matching.Rd
3a9546ce6447f5b6585b4bd58dc4c7f4 *man/as.phylo.Rd
219cff7b94e167e18f4f9cd99cc8efc3 *man/as.phylo.formula.Rd
@@ -166,12 +168,13 @@ ad514b163e70bfbc11dfd34a450799f8 *man/balance.Rd
868d03a447b2375bd9e99cde899cedea *man/base.freq.Rd
524a1163eac56d1fc7f65bcd6c74a8d0 *man/bd.ext.Rd
5f1e61abe5908708f503cb2626a1b5e6 *man/bd.time.Rd
+f929bc1b6391c57a6b0099c4561fd7be *man/binaryPGLMM.Rd
814a2437b237e083dfb91b404a8ede60 *man/bind.tree.Rd
822558d4f7ee04257b2c903c53ec4344 *man/bionj.Rd
71a008cfe65c4f524a5b66e68bbf81ab *man/bird.families.Rd
0e41770e1e6d0b8d90c4cf51049213cb *man/bird.orders.Rd
ef1c15d5d93410c21179997431112209 *man/birthdeath.Rd
-7569c1d30cd258f7ce414e2f42e7bc7a *man/boot.phylo.Rd
+428a6f9e08291416b6b96d945d734aac *man/boot.phylo.Rd
5a64b90d3a6c7a8204946b00f45f4cfc *man/branching.times.Rd
7dad9d632d914de8122f2455f0762546 *man/c.phylo.Rd
9d3f9614732671a610cc8a37454885e2 *man/carnivora.Rd
@@ -181,9 +184,9 @@ b40dd77fd2529626e983ee0300be6781 *man/chronoMPL.Rd
c1f01c6200b2f1e2901d45d40daae404 *man/chronopl.Rd
506d0332fb092ab87ca7674faef63ab7 *man/chronos.Rd
cb0cdda41058b55d02b87df14c86f85e *man/clustal.Rd
-479f0a9817a0cfc257badeca23660ef0 *man/coalescent.intervals.Rd
+866af6e8d769b3d6972ef8e1ac849a12 *man/coalescent.intervals.Rd
0841b19345d6f61e97852f6c810cccfd *man/collapse.singles.Rd
-08dfcbd8455de6b872a63883e38f5813 *man/collapsed.intervals.Rd
+bff5a7826f5a39767601e32ceb776247 *man/collapsed.intervals.Rd
301f271dc131de2efc3294d31f03afed *man/compar.cheverud.Rd
4d8ee141d7b6b323ef5ee9446000ae32 *man/compar.gee.Rd
4317f601ba2eef4e5730dc4274172f0f *man/compar.lynch.Rd
@@ -199,15 +202,16 @@ a210fe55aacb847936004458d4089b6d *man/corBrownian.Rd
673c9738c6ff3e237e11d35d62629bcc *man/corGrafen.Rd
628a4d6339f82e7ef3debe96d19ffa02 *man/corMartins.Rd
202f26f668a8fc7973d50696dcdd6cd4 *man/corPagel.Rd
+e259ee771509883d3afe8eafd41e9cb1 *man/corphylo.Rd
d7eb9b4fdf7036e82b5964bfd85e5e36 *man/correlogram.formula.Rd
c199605f9d353b303acad4896f9b39a5 *man/cynipids.Rd
34164e368efd0d5d961fe62e9ede75e8 *man/dbd.Rd
-e7703bd85471f644eba7aad66a1c5ad8 *man/def.Rd
-e5fd45ae77e515b42d180a8c5369987d *man/del.gaps.Rd
+006b1175da065de114ad006fe1e581d5 *man/def.Rd
+f4e8a12b122ff902a1bce3ec49e17da6 *man/del.gaps.Rd
fbcd1d4bcf74e21fc93c195c7af3db98 *man/delta.plot.Rd
e9a4068c29af95a6c2474192b3e1af1b *man/dist.dna.Rd
38011e81d28a120d88eead09e62c154a *man/dist.gene.Rd
-a93e4fe0f50b5b8ded1b47cb9c76ec1f *man/dist.topo.Rd
+e813d96ce8d458a535fffbb072701300 *man/dist.topo.Rd
c7cc398115be066740ca4fb037394727 *man/diversi.gof.Rd
d646ea0343999bd0e38e86dcf6c12018 *man/diversi.time.Rd
0003bd93b9418ba7a7a4f89faf2f7338 *man/diversity.contrast.test.Rd
@@ -221,13 +225,13 @@ eea313e8ee32597b4cec120d23113642 *man/gammaStat.Rd
18012643a904f657fc5f5896f3d14054 *man/howmanytrees.Rd
86c49d080fdffd614d8056021e91cc55 *man/identify.phylo.Rd
46c7c675e02d355a25d8c1c02deffcd7 *man/image.DNAbin.Rd
-4bb168908d445690105ad4c00661b93b *man/is.binary.tree.Rd
+9957425ef61231ba72538e771c6565b5 *man/is.binary.tree.Rd
80a5a228a43679edf75903590253a875 *man/is.compatible.Rd
d2de8fd9549ef01a1dddeb726dd77fcf *man/is.monophyletic.Rd
9d09e72a00da6b7b3f967f76166077cd *man/is.ultrametric.Rd
3f6ff8340b6c9770a1f4d2fed72a295d *man/kronoviz.Rd
2afa6305e48e4c47a7d94d46401f85a3 *man/ladderize.Rd
-30290ed7800e405b5e82e23af384c6af *man/landplants.Rd
+23db31affc0dc60f6c772e850bb640c3 *man/landplants.Rd
4b86460e4e5d993abf13c99b9744dbd6 *man/lmorigin.Rd
21bb31db5bcc8d8650467ef73fe0c0d3 *man/ltt.plot.Rd
3d81c146cea3555ac60a7dede5fbb882 *man/makeLabel.Rd
@@ -247,7 +251,7 @@ dbcdc231c2f1f221818307bb33486159 *man/mrca.Rd
00fb7ade93c2dd887be27e3dad8e2115 *man/mvr.Rd
3df9e16b8a09df3f1dba5c4327a635fc *man/nj.Rd
9ea7d5899a190171de4165e3240de62e *man/njs.Rd
-8f30d041ced38534477c2fac72fb1abf *man/node.depth.Rd
+a589b4cc04505185dc9ef1817c7ae304 *man/node.depth.Rd
ba754eba67181a1c39d417ffd9d23b76 *man/nodelabels.Rd
447ae03684ff56a4a24932aec182acf3 *man/nodepath.Rd
4f8b4cf2a725a8001e49390056181bef *man/opsin.Rd
@@ -266,21 +270,21 @@ b24438c42cea969302ec6ba61002426e *man/print.phylo.Rd
81f756fdf2ec4c968095595dece8338f *man/rTraitMult.Rd
71b55f1d67b63610555f316ac454fb58 *man/read.GenBank.Rd
9e31bccfa08c1d2da2180cd8f50ab250 *man/read.caic.Rd
-37d56f1aa48a2d9469d2e07d6dedc4d5 *man/read.dna.Rd
+0c91321b05d51c501e0f5c51de78538c *man/read.dna.Rd
5c9aed0b6f0a7e0f0fac9b5a59ac732f *man/read.nexus.Rd
bc02e36c51d67074e661468993ed359b *man/read.nexus.data.Rd
-b1af1fd1362b66f6eae928bef96b38e8 *man/read.tree.Rd
+3c8b53f5af4fdd36e7c0b449412d5a77 *man/read.tree.Rd
48d08588177df45a1e53dede1dc56bae *man/reconstruct.Rd
61e0bc86ae5a1f04ee6e88cdc92a2089 *man/reorder.phylo.Rd
af19f262b468dfd7c12448d9d8085331 *man/richness.yule.test.Rd
-a3cf7db0d79645f61aa15634b21a5c92 *man/rlineage.Rd
-aa7d8a43e8456fe53480c023b7ef1648 *man/root.Rd
+d919cbc29acafcefaa9a5b2594da6a8c *man/rlineage.Rd
+bb3a58966e6d9d72a4b885bbe9fb935f *man/root.Rd
6b97ea2fd96bf948fca764deab3a6e76 *man/rotate.Rd
1a7314ee0aa76712737e76162df6cc74 *man/rtree.Rd
e53b4b3f7e009289c628e16143b3e8b4 *man/rtt.Rd
-2b2faabf81e29f320eff5adde7e6d75d *man/seg.sites.Rd
-ce6f293b22e9b0fd36f38c67b5f40e6e *man/skyline.Rd
-74d4e9d6fc374c26807825795dff1e0e *man/skylineplot.Rd
+5010f7e12af0865bd3e4112349d601db *man/seg.sites.Rd
+f125fe172ee83b8bb4adaede4b4d3b43 *man/skyline.Rd
+bf851aa61b6310afa2ae1358c189dad7 *man/skylineplot.Rd
0cede7cdef45fc9a8b993f938a32ce67 *man/slowinskiguyer.test.Rd
a0db170675b1c265f1f4e7d56e88d7cb *man/speciesTree.Rd
8fc51ff26614e7c713be2570d9af47b6 *man/stree.Rd
@@ -296,10 +300,10 @@ dbb269e680caf8c722cf53c9b5ce7ace *man/triangMtd.Rd
5c459720196654a10cfff0390bffc14f *man/vcv.phylo.Rd
f7ce7760e913c10f1cb7be05315b39fc *man/vcv2phylo.Rd
a5d3cbf19df84d0f1e4e3a4650945cbf *man/weight.taxo.Rd
-19d511057993a6335549d8e84f6ef86a *man/where.Rd
+5d20c7995514d810a36e6ecf4c5faa87 *man/where.Rd
ef9658f5343bdcbfc3673c7f936934f5 *man/which.edge.Rd
-f569d0058bc598bb92c42ef99a58f272 *man/woodmouse.Rd
-510a3ad17baa847a3f581e9d1babb534 *man/write.dna.Rd
+9b83a148295dcfcb3994fdb6eaab0174 *man/woodmouse.Rd
+f613518d0ecd88b4a80c779b0453ad7f *man/write.dna.Rd
b9fe1836be61c56a521386b00021fa67 *man/write.nexus.Rd
80e6bcbf556ad15770c7157aef3d28d8 *man/write.nexus.data.Rd
b2302d65196726ab4cc573a22ee2d8e3 *man/write.tree.Rd
@@ -320,7 +324,7 @@ d343190ba59ef1feac9c96752f210687 *src/ape.c
f782a97466ed8008d31afd7a3e9d4ea9 *src/bipartition.c
a5d4f692aca36f6b50e123b443792fa1 *src/bitsplits.c
81e4ad60b418966cdd09e608e074bb28 *src/delta_plot.c
-db663ebe2c21f6824447223b7322de59 *src/dist_dna.c
+00d3d2104b1b606b17d4cc636129d2e7 *src/dist_dna.c
2a6a59a3a220eb855a3c48bc73dc54ba *src/dist_nodes.c
005ab69356525fbbc1b69950341308c2 *src/ewLasso.c
afa27403972145d083d7f7f8b2536b98 *src/heap.c
diff --git a/NAMESPACE b/NAMESPACE
index f2627e1..ded7b7f 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -21,14 +21,19 @@ S3method(deviance, ace)
S3method(logLik, ace)
S3method(print, ace)
+S3method(print, binaryPGLMM)
+
S3method(as.prop.part, bitsplits)
S3method(is.compatible, bitsplits)
S3method(print, bitsplits)
+S3method(sort, bitsplits)
S3method(drop1, compar.gee)
S3method(predict, compar.gee)
S3method(print, compar.gee)
+S3method(print, corphylo)
+
S3method("[", DNAbin)
S3method(as.character, DNAbin)
S3method(as.list, DNAbin)
diff --git a/R/CDF.birth.death.R b/R/CDF.birth.death.R
index 8dadfdb..c632528 100644
--- a/R/CDF.birth.death.R
+++ b/R/CDF.birth.death.R
@@ -1,9 +1,8 @@
-## CDF.birth.death.R (2014-03-03)
+## CDF.birth.death.R (2015-04-10)
-## Functions to Simulate and Fit
-## Time-Dependent Birth-Death Models
+## Functions to Simulate and Fit Time-Dependent Birth-Death Models
-## Copyright 2010-2014 Emmanuel Paradis
+## Copyright 2010-2015 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
@@ -232,7 +231,7 @@ rlineage <-
},{ # case 7:
Foo <- function(t, x) {
if (t == x) return(0)
- 1 - exp(-(birth*(x - t) + DEATH(x) - DEATH(t) ))
+ 1 - exp(-(birth*(x - t) + DEATH(x) - DEATH(t)))
}
})
@@ -517,3 +516,168 @@ LTT <- function(birth = 0.1, death = 0, N = 100, Tmax = 50, PI = 95,
lines(c(Time, NA, Time), c(Flow, NA, Fup),
col = pi.style[[1]], lwd = pi.style[[2]], lty = pi.style[[3]])
}
+
+rphylo <-
+ function(n, birth, death, BIRTH = NULL, DEATH = NULL,
+ T0 = 50, fossils = FALSE, eps = 1e-6)
+{
+ case <- .getCase(birth, death, BIRTH, DEATH)
+
+ ## Foo(): CDF of the times to event (speciation or extinction)
+ ## rTimeToEvent(): generate a random time to event by the inverse method
+
+ switch(case, { # case 1:
+ rTimeToEvent <- function(t) t - rexp(1, N * (birth + death)) # much faster than using Foo()
+ speORext <- function(t) birth/(birth + death)
+ ## Foo <- function(t, x)
+ ## 1 - exp(-N*(birth + death)*(x - t))
+ },{ # case 2:
+ Foo <- function(t, x) {
+ if (t == x) return(0)
+ 1 - exp(-integrate(function(t) birth(t) + death(t),
+ t, x)$value * N)
+ }
+ speORext <- function(t) birth(t)/(birth(t) + death(t))
+ },{ # case 3:
+ Foo <- function(t, x) {
+ if (t == x) return(0)
+ 1 - exp(-N*(BIRTH(x) - BIRTH(t) + DEATH(x) - DEATH(t)))
+ }
+ speORext <- function(t) birth(t)/(birth(t) + death(t))
+ },{ # case 4:
+ Foo <- function(t, x) {
+ if (t == x) return(0)
+ 1 - exp(-N*(integrate(function(t) birth(t), t, x)$value
+ + death*(x - t)))
+ }
+ speORext <- function(t) birth(t)/(birth(t) + death)
+ },{ # case 5:
+ Foo <- function(t, x) {
+ if (t == x) return(0)
+ 1 - exp(-N*(birth*(x - t) +
+ integrate(function(t) death(t), t, x)$value))
+ }
+ speORext <- function(t) birth/(birth + death(t))
+ },{ # case 6:
+ Foo <- function(t, x) {
+ if (t == x) return(0)
+ 1 - exp(-N*(BIRTH(x) - BIRTH(t) + death*(x - t)))
+ }
+ speORext <- function(t) birth(t)/(birth(t) + death)
+ },{ # case 7:
+ Foo <- function(t, x) {
+ if (t == x) return(0)
+ 1 - exp(-N*(birth*(x - t) + DEATH(x) - DEATH(t)))
+ }
+ speORext <- function(t) birth/(birth + death(t))
+ })
+
+ if (case != 1) {
+ rTimeToEvent <- function(t)
+ {
+ P <- runif(1)
+ inc <- 10
+ x <- t - inc
+ while (inc > eps) {
+ if (Foo(x, t) > P) {
+ x <- x + inc
+ inc <- inc/10
+ }
+ x <- x - inc
+ }
+ x
+ }
+ }
+
+ storage.mode(n) <- "integer"
+ N <- n
+ t <- T0
+ j <- 0L # number of edges already created
+ POOL <- seq_len(N) # initial pool (only tips at start)
+
+if (!fossils) {
+ Nedge <- 2L * N - 2L
+ nextnode <- 2L * N - 1L
+ e1 <- integer(Nedge)
+ e2 <- integer(Nedge)
+ TIME <- numeric(nextnode) # record the times
+ TIME[POOL] <- T0 # the times of the n tips are the present time
+
+ while (j < Nedge) {
+ X <- rTimeToEvent(t)
+ ## is the event a speciation or an extinction?
+ Y <- rbinom(1, size = 1, prob = speORext(X))
+ if (Y) { # speciation
+ i <- sample.int(N, 2)
+ fossil <- POOL[i] == 0
+ if (any(fossil)) {
+ ## we drop the fossil lineage, or the first one if both are fossils
+ POOL <- POOL[-i[which(fossil)[1]]]
+ } else { # create a node and an edge
+ j <- j + 2L
+ k <- c(j - 1, j)
+ e1[k] <- nextnode
+ e2[k] <- POOL[i]
+ TIME[nextnode] <- X
+ POOL <- c(POOL[-i], nextnode)
+ nextnode <- nextnode - 1L
+ }
+ N <- N - 1L
+ } else { # extinction => create a tip, store it in POOL but don't create an edge
+ ## fossil lineages are numbered 0 to find them if Y = 1
+ N <- N + 1L
+ POOL <- c(POOL, 0L)
+ }
+ t <- X
+ }
+
+ Nnode <- n - 1L
+
+} else { # fossils = TRUE
+ nextnode <- -1L # nodes are numbered with negatives
+ nexttip <- N + 1L # tips are numbered with positives
+ e1 <- integer(1e5)
+ e2 <- integer(1e5)
+ time.tips <- numeric(1e5) # accessed with positive indices
+ time.nodes <- numeric(1e5) # accessed with negative indices
+ time.tips[POOL] <- T0 # the times of the n living tips are the present time
+
+ while (N > 1) {
+ X <- rTimeToEvent(t)
+ ## is the event a speciation or an extinction?
+ Y <- rbinom(1, size = 1, prob = speORext(X))
+ if (Y) { # speciation => create a node
+ i <- sample.int(N, 2)
+ j <- j + 2L
+ k <- c(j - 1, j)
+ e1[k] <- nextnode
+ e2[k] <- POOL[i]
+ time.nodes[-nextnode] <- X
+ POOL <- c(POOL[-i], nextnode)
+ nextnode <- nextnode - 1L
+ N <- N - 1L
+ } else { # extinction => create a tip
+ N <- N + 1L
+ time.tips[nexttip] <- X
+ POOL <- c(POOL, nexttip)
+ nexttip <- nexttip + 1L
+ }
+ t <- X
+ }
+
+ n <- nexttip - 1L # update n
+ Nnode <- n - 1L
+ EDGE <- seq_len(j)
+ e1 <- e1[EDGE]
+ e2 <- e2[EDGE]
+ e1 <- e1 + n + Nnode + 1L # e1 has only nodes...
+ NODES <- e2 < 0 # ... so this is needed only on e2
+ e2[NODES] <- e2[NODES] + n + Nnode + 1L
+ ## concatenate the vectors of times after dropping the extra 0's:
+ TIME <- c(time.tips[seq_len(n)], rev(time.nodes[seq_len(Nnode)]))
+}
+ structure(list(edge = cbind(e1, e2, deparse.level = 0),
+ edge.length = TIME[e2] - TIME[e1],
+ tip.label = paste0("t", seq_len(n)),
+ Nnode = Nnode), class = "phylo")
+}
diff --git a/R/DNA.R b/R/DNA.R
index 4fed6ce..8306eb4 100644
--- a/R/DNA.R
+++ b/R/DNA.R
@@ -1,8 +1,8 @@
-## DNA.R (2014-08-05)
+## DNA.R (2015-05-28)
## Manipulations and Comparisons of DNA Sequences
-## Copyright 2002-2014 Emmanuel Paradis
+## Copyright 2002-2015 Emmanuel Paradis, 2015 Klaus Schliep
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
@@ -36,6 +36,17 @@ del.gaps <- function(x)
x
}
+del.colgapsonly <- function(x)
+{
+ if (!inherits(x, "DNAbin")) x <- as.DNAbin(x)
+ if (!is.matrix(x))
+ stop("DNA sequences not in a matrix")
+ foo <- function(x) all(x == 4)
+ i <- which(apply(x, 2, foo))
+ if (length(i)) x <- x[, -i]
+ x
+}
+
as.alignment <- function(x)
{
if (is.list(x)) n <- length(x)
@@ -231,14 +242,10 @@ as.DNAbin <- function(x, ...) UseMethod("as.DNAbin")
._bs_ <- c(136, 72, 40, 24, 192, 160, 144, 96, 80,
48, 224, 176, 208, 112, 240, 4, 2)
+## by Klaus:
as.DNAbin.character <- function(x, ...)
{
- n <- length(x)
- ans <- raw(n)
- for (i in 1:15)
- ans[which(x == ._cs_[i])] <- as.raw(._bs_[i])
- ans[which(x == "-")] <- as.raw(4)
- ans[which(x == "?")] <- as.raw(2)
+ ans <- as.raw(._bs_)[match(x, ._cs_)]
if (is.matrix(x)) {
dim(ans) <- dim(x)
dimnames(ans) <- dimnames(x)
@@ -268,11 +275,7 @@ as.DNAbin.list <- function(x, ...)
as.character.DNAbin <- function(x, ...)
{
f <- function(xx) {
- ans <- character(length(xx))
- for (i in 1:15)
- ans[which(xx == ._bs_[i])] <- ._cs_[i]
- ans[which(xx == 4)] <- "-"
- ans[which(xx == 2)] <- "?"
+ ans <- ._cs_[match(as.numeric(xx), ._bs_)]
if (is.matrix(xx)) {
dim(ans) <- dim(xx)
dimnames(ans) <- dimnames(xx)
diff --git a/R/PGLS.R b/R/PGLS.R
index ed7bc2b..c7dba5c 100644
--- a/R/PGLS.R
+++ b/R/PGLS.R
@@ -1,8 +1,8 @@
-## PGLS.R (2012-02-09)
+## PGLS.R (2014-12-11)
## Phylogenetic Generalized Least Squares
-## Copyright 2004 Julien Dutheil, and 2006-2012 Emmanuel Paradis
+## Copyright 2004 Julien Dutheil, and 2006-2014 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
@@ -92,7 +92,6 @@ corMatrix.corBrownian <-
stop("object has not been initialized.")
tree <- attr(object, "tree")
mat <- vcv.phylo(tree, corr = corr)
- n <- dim(mat)[1]
## reorder matrix:
index <- attr(object, "index")
mat[index, index]
@@ -109,7 +108,6 @@ corMatrix.corMartins <-
dist <- cophenetic.phylo(tree)
mat <- exp(-object[1] * dist)
if (corr) mat <- cov2cor(mat)
- n <- dim(mat)[1]
## reorder matrix:
index <- attr(object, "index")
mat[index, index]
@@ -125,7 +123,6 @@ corMatrix.corGrafen <-
tree <- compute.brlen(attr(object, "tree"),
method = "Grafen", power = exp(object[1]))
mat <- vcv.phylo(tree, corr = corr)
- n <- dim(mat)[1]
## reorder matrix:
index <- attr(object, "index")
mat[index, index]
diff --git a/R/ace.R b/R/ace.R
index fd19340..54a2be6 100644
--- a/R/ace.R
+++ b/R/ace.R
@@ -1,8 +1,8 @@
-## ace.R (2014-05-27)
+## ace.R (2015-05-01)
## Ancestral Character Estimation
-## Copyright 2005-2014 Emmanuel Paradis and Ben Bolker
+## Copyright 2005-2015 Emmanuel Paradis and Ben Bolker
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
@@ -227,10 +227,12 @@ ace <-
if (is.na(dev)) Inf else dev
}
} else {
- E <- if (use.expm) {
- library(expm)
- get("expm", "package:expm") # to avoid Matrix::expm
- } else matexpo
+ if (!requireNamespace("expm", quietly = TRUE) && use.expm) {
+ warning("package 'expm' not available; using function 'matexpo' from 'ape'")
+ use.expm <- FALSE
+ }
+ E <- if (use.expm) expm::expm # to avoid Matrix::expm
+ else matexpo
dev <- function(p, output.liks = FALSE) {
if (any(is.nan(p)) || any(is.infinite(p))) return(1e50)
comp <- numeric(nb.tip + nb.node) # from Rich FitzJohn
diff --git a/R/as.bitsplits.R b/R/as.bitsplits.R
index 7201bd8..8dcca9f 100644
--- a/R/as.bitsplits.R
+++ b/R/as.bitsplits.R
@@ -1,4 +1,4 @@
-## as.bitsplits.R (2014-01-02)
+## as.bitsplits.R (2014-12-11)
## Conversion Among Split Classes
@@ -57,6 +57,14 @@ print.bitsplits <- function(x, ...)
cat(' ', length(x$freq), 'partitions\n')
}
+sort.bitsplits <- function(x, decreasing = FALSE, ...)
+{
+ o <- order(x$freq, decreasing = decreasing)
+ x$matsplit <- x$matsplit[, o]
+ x$freq <- x$freq[o]
+ x
+}
+
as.prop.part <- function(x) UseMethod("as.prop.part")
as.prop.part.bitsplits <- function(x)
diff --git a/R/as.phylo.R b/R/as.phylo.R
index b56e5c2..fa2f373 100644
--- a/R/as.phylo.R
+++ b/R/as.phylo.R
@@ -1,8 +1,8 @@
-## as.phylo.R (2013-08-13)
+## as.phylo.R (2015-02-06)
## Conversion Among Tree Objects
-## Copyright 2005-2013 Emmanuel Paradis
+## Copyright 2005-2015 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
@@ -29,8 +29,7 @@ new2old.phylo <- function(phy)
as.phylo <- function (x, ...)
{
- if (length(class(x)) == 1 && class(x) == "phylo")
- return(x)
+ if (identical(class(x), "phylo")) return(x)
UseMethod("as.phylo")
}
diff --git a/R/binaryPGLMM.R b/R/binaryPGLMM.R
new file mode 100644
index 0000000..98d834a
--- /dev/null
+++ b/R/binaryPGLMM.R
@@ -0,0 +1,293 @@
+## binaryPGLMM.R (2015-03-04)
+
+## Phylogenetic Generalized Linear Mixed Model for Binary Data
+
+## Copyright 2015 Anthony R. Ives
+
+## This file is part of the R-package `ape'.
+## See the file ../COPYING for licensing issues.
+
+binaryPGLMM <- function(formula, data = list(), phy, s2.init = 0.1, B.init = NULL,
+ tol.pql = 10^-6, maxit.pql = 200, maxit.reml = 100) {
+
+ # Begin pglmm.reml
+ pglmm.reml <- function(par, tinvW, tH, tVphy, tX) {
+ n <- dim(tX)[1]
+ p <- dim(tX)[2]
+
+ ss2 <- abs(Re(par))
+ Cd <- ss2 * tVphy
+ V <- tinvW + Cd
+
+ LL <- 10^10
+ if (sum(is.infinite(V)) == 0) { # & rcond(V) < 10^10) {
+ if (all(eigen(V)$values > 0)) { #if(rcond(V) > 10^-10 & all(eigen(V)$values > 0)) {
+ invV <- solve(V)
+ logdetV <- determinant(V)$modulus[1]
+ if (is.infinite(logdetV)) {
+ cholV <- chol(V)
+ logdetV <- 2 * sum(log(diag(chol(V))))
+ }
+ LL <- logdetV + t(tH) %*% invV %*% tH + determinant(t(tX) %*%
+ invV %*% tX)$modulus[1]
+ }
+ }
+ return(LL)
+ }
+ # End pglmm.reml
+
+ if (!inherits(phy, "phylo"))
+ stop("Object \"phy\" is not of class \"phylo\".")
+ if (is.null(phy$edge.length))
+ stop("The tree has no branch lengths.")
+ if (is.null(phy$tip.label))
+ stop("The tree has no tip labels.")
+ phy <- reorder(phy, "postorder")
+ n <- length(phy$tip.label)
+
+ mf <- model.frame(formula = formula, data = data)
+ if (nrow(mf) != length(phy$tip.label))
+ stop("Number of rows of the design matrix does not match with length of the tree.")
+ if (is.null(rownames(mf))) {
+ warning("No tip labels, order assumed to be the same as in the tree.\n")
+ data.names = phy$tip.label
+ } else data.names = rownames(mf)
+ order <- match(data.names, phy$tip.label)
+ if (sum(is.na(order)) > 0) {
+ warning("Data names do not match with the tip labels.\n")
+ rownames(mf) <- data.names
+ } else {
+ tmp <- mf
+ rownames(mf) <- phy$tip.label
+ mf[order, ] <- tmp[1:nrow(tmp), ]
+ }
+
+ X <- model.matrix(attr(mf, "terms"), data = mf)
+ y <- model.response(mf)
+ if (sum(!(y %in% c(0, 1)))) {
+ stop("PGLMM.binary requires a binary response (dependent variable).")
+ }
+ if (var(y) == 0) {
+ stop("The response (dependent variable) is always 0 or always 1.")
+ }
+
+ p <- ncol(X)
+ Vphy <- vcv(phy)
+ Vphy <- Vphy/max(Vphy)
+ Vphy/exp(determinant(Vphy)$modulus[1]/n)
+
+ # Compute initial estimates if not provided assuming no phylogeny
+ if (!is.null(B.init) & length(B.init) != p) {
+ warning("B.init not correct length, so computed B.init using glm()")
+ }
+ if (is.null(B.init) | (!is.null(B.init) & length(B.init) != p)) {
+ B.init <- t(matrix(glm(formula = formula, data = data, family = "binomial")$coefficients, ncol = p))
+ }
+ B <- B.init
+ s2 <- s2.init
+ b <- matrix(0, nrow = n)
+ beta <- rbind(B, b)
+ mu <- exp(X %*% B)/(1 + exp(X %*% B))
+
+ XX <- cbind(X, diag(1, nrow = n, ncol = n))
+ C <- s2 * Vphy
+
+ est.s2 <- s2
+ est.B <- B
+ oldest.s2 <- 10^6
+ oldest.B <- matrix(10^6, nrow = length(est.B))
+
+ iteration <- 0
+ exitflag <- 0
+ rcondflag <- 0
+ while (((t(est.s2 - oldest.s2) %*% (est.s2 - oldest.s2) > tol.pql^2) |
+ (t(est.B - oldest.B) %*% (est.B - oldest.B)/length(B) > tol.pql^2)) &
+ (iteration <= maxit.pql)) {
+
+ iteration <- iteration + 1
+ oldest.s2 <- est.s2
+ oldest.B <- est.B
+
+ est.B.m <- B
+ oldest.B.m <- matrix(10^6, nrow = length(est.B))
+ iteration.m <- 0
+
+ # mean component
+ while ((t(est.B.m - oldest.B.m) %*% (est.B.m - oldest.B.m)/length(B) >
+ tol.pql^2) & (iteration.m <= maxit.pql)) {
+
+ iteration.m <- iteration.m + 1
+ oldest.B.m <- est.B.m
+ invW <- diag(as.vector((mu * (1 - mu))^-1))
+ V <- invW + C
+
+ # This flags cases in which V has a very high condition number, which will cause solve() to fail.
+ if (sum(is.infinite(V)) > 0 | rcond(V) < 10^-10) {
+ rcondflag <- rcondflag + 1
+ B <- 0 * B.init + 0.001
+ b <- matrix(0, nrow = n)
+ beta <- rbind(B, b)
+ mu <- exp(X %*% B)/(1 + exp(X %*% B))
+ oldest.B.m <- matrix(10^6, nrow = length(est.B))
+ invW <- diag(as.vector((mu * (1 - mu))^-1))
+ V <- invW + C
+ }
+
+ invV <- solve(V)
+ Z <- X %*% B + b + (y - mu)/(mu * (1 - mu))
+ denom <- t(X) %*% invV %*% X
+ num <- t(X) %*% invV %*% Z
+ B <- as.matrix(solve(denom, num))
+
+ b <- C %*% invV %*% (Z - X %*% B)
+ beta <- rbind(B, b)
+ mu <- exp(XX %*% beta)/(1 + exp(XX %*% beta))
+
+ est.B.m <- B
+ }
+
+ # variance component
+ H <- Z - X %*% B
+ opt <- optim(fn = pglmm.reml, par = s2, tinvW = invW, tH = H, tVphy = Vphy,
+ tX = X, method = "BFGS", control = list(factr = 1e+12, maxit = maxit.reml))
+ s2 <- abs(opt$par)
+ C <- s2 * Vphy
+
+ est.s2 <- s2
+ est.B <- B
+ }
+ convergeflag <- "converged"
+ if (iteration >= maxit.pql | rcondflag >= 3) {
+ convergeflag <- "Did not converge; try increasing maxit.pql or starting with B.init values of .001"
+ }
+
+ converge.test.s2 <- (t(est.s2 - oldest.s2) %*% (est.s2 - oldest.s2))^0.5
+ converge.test.B <- (t(est.B - oldest.B) %*% (est.B - oldest.B))^0.5/length(est.B)
+
+ # Extract parameters
+ invW <- diag(as.vector((mu * (1 - mu))^-1))
+ V <- invW + C
+ invV <- solve(V)
+ Z <- X %*% B + b + (y - mu)/(mu * (1 - mu))
+ denom <- t(X) %*% invV %*% X
+ num <- t(X) %*% invV %*% Z
+ B <- solve(denom, num)
+ b <- C %*% invV %*% (Z - X %*% B)
+ beta <- rbind(B, b)
+ mu <- exp(XX %*% beta)/(1 + exp(XX %*% beta))
+ H <- Z - X %*% B
+
+ B.cov <- solve(t(X) %*% invV %*% X)
+ B.se <- as.matrix(diag(B.cov))^0.5
+ B.zscore <- B/B.se
+ B.pvalue <- 2 * pnorm(abs(B/B.se), lower.tail = FALSE)
+
+ LL <- opt$value
+ lnlike.cond.reml <- -0.5 * (n - p) * log(2 * pi) + 0.5 * determinant(t(X) %*%
+ X)$modulus[1] - 0.5 * LL
+
+ LL0 <- pglmm.reml(par = 0, tinvW = invW, tH = H, tVphy = Vphy, tX = X)
+ lnlike.cond.reml0 <- -0.5 * (n - p) * log(2 * pi) + 0.5 * determinant(t(X) %*%
+ X)$modulus[1] - 0.5 * LL0
+
+ P.H0.s2 <- pchisq(2 * (lnlike.cond.reml - lnlike.cond.reml0), df = 1,
+ lower.tail = F)/2
+
+ results <- list(formula = formula, B = B, B.se = B.se, B.cov = B.cov,
+ B.zscore = B.zscore, B.pvalue = B.pvalue, s2 = s2, P.H0.s2 = P.H0.s2,
+ mu = mu, b = b, B.init = B.init, X = X, H = H, VCV = Vphy, V = V,
+ convergeflag = convergeflag, iteration = iteration, converge.test.s2 = converge.test.s2,
+ converge.test.B = converge.test.B, rcondflag = rcondflag)
+ class(results) <- "binaryPGLMM"
+ results
+}
+
+######################################################
+######################################################
+# binaryPGLMM.sim
+######################################################
+######################################################
+
+binaryPGLMM.sim <- function(formula, data = list(), phy, s2 = NULL, B = NULL,
+ nrep = 1) {
+
+ if (!inherits(phy, "phylo"))
+ stop("Object \"phy\" is not of class \"phylo\".")
+ if (is.null(phy$edge.length))
+ stop("The tree has no branch lengths.")
+ if (is.null(phy$tip.label))
+ stop("The tree has no tip labels.")
+ phy <- reorder(phy, "postorder")
+ n <- length(phy$tip.label)
+
+ mf <- model.frame(formula = formula, data = data)
+ if (nrow(mf) != length(phy$tip.label))
+ stop("Number of rows of the design matrix does not match with length of the tree.")
+ if (is.null(rownames(mf))) {
+ warning("No tip labels, order assumed to be the same as in the tree.\n")
+ data.names = phy$tip.label
+ } else data.names = rownames(mf)
+ order <- match(data.names, phy$tip.label)
+ if (sum(is.na(order)) > 0) {
+ warning("Data names do not match with the tip labels.\n")
+ rownames(mf) <- data.names
+ } else {
+ tmp <- mf
+ rownames(mf) <- phy$tip.label
+ mf[order, ] <- tmp[1:nrow(tmp), ]
+ }
+
+ if (is.null(s2))
+ stop("You must specify s2")
+ if (is.null(B))
+ stop("You must specify B")
+ X <- model.matrix(attr(mf, "terms"), data = mf)
+
+ n <- nrow(X)
+ p <- ncol(X)
+ V <- vcv(phy)
+ V <- V/max(V)
+ V <- vcv(phy)
+ V <- V/max(V)
+ V/exp(determinant(V)$modulus[1]/n)
+ V <- s2 * V
+
+ if (s2 > 10^-8) {
+ iD <- t(chol(V))
+ } else {
+ iD <- matrix(0, nrow = n, ncol = n)
+ }
+ Y <- matrix(0, nrow = n, ncol = nrep)
+ y <- matrix(0, nrow = n, ncol = nrep)
+ for (i in 1:nrep) {
+ y[, i] <- X %*% B + iD %*% rnorm(n = n)
+ p <- 1/(1 + exp(-y[, i]))
+ Y[, i] <- as.numeric(runif(n = n) < p)
+ }
+ results <- list(Y = Y, y = y, X = X, s2 = s2, B = B, V = V)
+ return(results)
+}
+
+######################################################
+######################################################
+# print.binaryPGLMM
+######################################################
+######################################################
+
+print.binaryPGLMM <- function(x, digits = max(3, getOption("digits") - 3),
+ ...) {
+
+ cat("\n\nCall:")
+ print(x$formula)
+ cat("\n")
+
+ cat("Random effect (phylogenetic signal s2):\n")
+ w <- data.frame(s2 = x$s2, Pr = x$P.H0.s2)
+ print(w, digits = digits)
+
+ cat("\nFixed effects:\n")
+ coef <- data.frame(Value = x$B, Std.Error = x$B.se, Zscore = x$B.zscore,
+ Pvalue = x$B.pvalue)
+ printCoefmat(coef, P.values = TRUE, has.Pvalue = TRUE)
+ cat("\n")
+}
diff --git a/R/compar.gee.R b/R/compar.gee.R
index 5edec9b..a996750 100644
--- a/R/compar.gee.R
+++ b/R/compar.gee.R
@@ -1,8 +1,8 @@
-## compar.gee.R (2013-12-19)
+## compar.gee.R (2015-05-01)
## Comparative Analysis with GEEs
-## Copyright 2002-2013 Emmanuel Paradis
+## Copyright 2002-2015 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
@@ -11,7 +11,8 @@ compar.gee <-
function(formula, data = NULL, family = gaussian, phy,
corStruct, scale.fix = FALSE, scale.value = 1)
{
- library(gee)
+ if (requireNamespace("gee", quietly = TRUE)) gee <- gee::gee
+ else stop("package 'gee' not available")
if (!missing(corStruct)) {
if (!missing(phy))
diff --git a/R/corphylo.R b/R/corphylo.R
new file mode 100644
index 0000000..9ec1df4
--- /dev/null
+++ b/R/corphylo.R
@@ -0,0 +1,355 @@
+## corphylo.R (2015-05-01)
+
+## Ancestral Character Estimation
+
+## Copyright 2015 Anthony R. Ives
+
+## This file is part of the R-package `ape'.
+## See the file ../COPYING for licensing issues.
+
+corphylo <- function(X, U = list(), SeM = NULL, phy = NULL, REML = TRUE, method = c("Nelder-Mead", "SANN"),
+ constrain.d = FALSE, reltol = 10^-6, maxit.NM = 1000, maxit.SA = 1000, temp.SA = 1, tmax.SA = 1, verbose = FALSE) {
+
+ # Begin corphylo.LL
+ corphylo.LL <- function(par, XX, UU, MM, tau, Vphy, REML, constrain.d, verbose) {
+
+ n <- nrow(X)
+ p <- ncol(X)
+
+ L.elements <- par[1:(p + p * (p - 1)/2)]
+ L <- matrix(0, nrow = p, ncol = p)
+ L[lower.tri(L, diag = T)] <- L.elements
+ R <- t(L) %*% L
+
+ if (constrain.d == TRUE) {
+ logit.d <- par[(p + p * (p - 1)/2 + 1):length(par)]
+ if (max(abs(logit.d)) > 10)
+ return(10^10)
+ d <- 1/(1 + exp(-logit.d))
+ } else {
+ d <- par[(p + p * (p - 1)/2 + 1):length(par)]
+ if (max(d) > 10)
+ return(10^10)
+ }
+
+ # OU transform
+ C <- matrix(0, nrow = p * n, ncol = p * n)
+ for (i in 1:p) for (j in 1:p) {
+ Cd <- (d[i]^tau * (d[j]^t(tau)) * (1 - (d[i] * d[j])^Vphy))/(1 - d[i] * d[j])
+ C[(n * (i - 1) + 1):(i * n), (n * (j - 1) + 1):(j * n)] <- R[i, j] * Cd
+ }
+
+ V <- C + diag(as.numeric(MM))
+ if (is.nan(rcond(V)) || rcond(V) < 10^-10)
+ return(10^10)
+ iV <- solve(V)
+ denom <- t(UU) %*% iV %*% UU
+ if (is.nan(rcond(denom)) || rcond(denom) < 10^-10)
+ return(10^10)
+ num <- t(UU) %*% iV %*% XX
+ B <- solve(denom, num)
+ B <- as.matrix(B)
+ H <- XX - UU %*% B
+
+ logdetV <- -determinant(iV)$modulus[1]
+ if (is.infinite(logdetV))
+ return(10^10)
+
+ if (REML == TRUE) {
+ # REML likelihood function
+ LL <- 0.5 * (logdetV + determinant(t(UU) %*% iV %*% UU)$modulus[1] + t(H) %*% iV %*% H)
+ } else {
+ # ML likelihood function
+ LL <- 0.5 * (logdetV + t(H) %*% iV %*% H)
+ }
+
+ if (verbose == T)
+ show(c(as.numeric(LL), par))
+ return(as.numeric(LL))
+ }
+ # End corphylo.LL
+
+ # Main program
+ if (!inherits(phy, "phylo"))
+ stop("Object \"phy\" is not of class \"phylo\".")
+ if (is.null(phy$edge.length))
+ stop("The tree has no branch lengths.")
+ if (is.null(phy$tip.label))
+ stop("The tree has no tip labels.")
+ phy <- reorder(phy, "postorder")
+ n <- length(phy$tip.label)
+
+ # Input X
+ if (dim(X)[1] != n)
+ stop("Number of rows of the data matrix does not match the length of the tree.")
+ if (is.null(rownames(X))) {
+ warning("No tip labels on X; order assumed to be the same as in the tree.\n")
+ data.names = phy$tip.label
+ } else data.names = rownames(X)
+ order <- match(data.names, phy$tip.label)
+ if (sum(is.na(order)) > 0) {
+ warning("Data names do not match with the tip labels.\n")
+ rownames(X) <- data.names
+ } else {
+ temp <- X
+ rownames(X) <- phy$tip.label
+ X[order, ] <- temp[1:nrow(temp), ]
+ }
+ p <- dim(X)[2]
+
+ # Input SeM
+ if (!is.null(SeM)) {
+ if (dim(SeM)[1] != n)
+ stop("Number of rows of the SeM matrix does not match the length of the tree.")
+ if (is.null(rownames(SeM))) {
+ warning("No tip labels on SeM; order assumed to be the same as in the tree.\n")
+ data.names = phy$tip.label
+ } else data.names = rownames(SeM)
+ order <- match(data.names, phy$tip.label)
+ if (sum(is.na(order)) > 0) {
+ warning("SeM names do not match with the tip labels.\n")
+ rownames(SeM) <- data.names
+ } else {
+ temp <- SeM
+ rownames(SeM) <- phy$tip.label
+ SeM[order, ] <- temp[1:nrow(temp), ]
+ }
+ } else {
+ SeM <- matrix(0, nrow = n, ncol = p)
+ }
+
+ # Input U
+ if (length(U) > 0) {
+ if (length(U) != p)
+ stop("Number of elements of list U does not match the number of columns in X.")
+
+ for (i in 1:p) {
+ if (!is.null(U[[i]])){
+ if (dim(U[[i]])[1] != n)
+ stop("Number of rows of an element of U does not match the tree.")
+ if (is.null(rownames(U[[i]]))) {
+ warning("No tip labels on U; order assumed to be the same as in the tree.\n")
+ data.names = phy$tip.label
+ } else data.names = rownames(U[[i]])
+ order <- match(data.names, phy$tip.label)
+ if (sum(is.na(order)) > 0) {
+ warning("U names do not match with the tip labels.\n")
+ rownames(U[[i]]) <- data.names
+ } else {
+ temp <- U[[i]]
+ rownames(U[[i]]) <- phy$tip.label
+ U[[i]][order, ] <- temp[1:nrow(temp), ]
+ }
+ } else {
+ U[[i]] <- matrix(0, nrow=n, ncol=1)
+ rownames(U[[i]]) <- phy$tip.label
+ }
+ }
+ }
+
+ # Standardize all variables
+ Xs <- X
+ for (i in 1:p) Xs[, i] <- (X[, i] - mean(X[, i]))/sd(X[, i])
+
+ if (!is.null(SeM)) {
+ SeMs <- SeM
+ for (i in 1:p) SeMs[, i] <- SeM[, i]/sd(X[, i])
+ }
+
+ if (length(U) > 0) {
+ Us <- U
+ for (i in 1:p) for (j in 1:ncol(U[[i]])) {
+ if (sd(U[[i]][, j]) > 0) {
+ Us[[i]][, j] <- (U[[i]][, j] - mean(U[[i]][, j]))/sd(U[[i]][, j])
+ } else {
+ Us[[i]][, j] <- U[[i]][, j] - mean(U[[i]][, j])
+ }
+ }
+ }
+
+ # Set up matrices
+ Vphy <- vcv(phy)
+ Vphy <- Vphy/max(Vphy)
+ Vphy <- Vphy/exp(determinant(Vphy)$modulus[1]/n)
+
+ XX <- matrix(as.matrix(Xs), ncol = 1)
+ MM <- matrix(as.matrix(SeMs^2), ncol = 1)
+
+ UU <- kronecker(diag(p), matrix(1, nrow = n, ncol = 1))
+ if (length(U) > 0) {
+ zeros <- 0 * (1:p)
+ for (i in 1:p) {
+ dd <- zeros
+ dd[i] <- 1
+ u <- kronecker(dd, as.matrix(Us[[i]]))
+ for (j in 1:dim(u)[2]) if (sd(u[, j]) > 0)
+ UU <- cbind(UU, u[, j])
+ }
+ }
+
+ # Compute initial estimates assuming no phylogeny if not provided
+ if (length(U) > 0) {
+ eps <- matrix(nrow = n, ncol = p)
+ for (i in 1:p) {
+ if (ncol(U[[i]]) > 0) {
+ u <- as.matrix(Us[[i]])
+ z <- lm(Xs[, i] ~ u)
+ eps[, i] <- resid(z)
+ } else {
+ eps[, i] <- Xs[, i] - mean(Xs[, i])
+ }
+ }
+ L <- t(chol(cov(eps)))
+ } else {
+ L <- t(chol(cov(Xs)))
+ }
+ L.elements <- L[lower.tri(L, diag = T)]
+ par <- c(L.elements, array(0.5, dim = c(1, p)))
+
+ tau <- matrix(1, nrow = n, ncol = 1) %*% diag(Vphy) - Vphy
+
+ if (method == "Nelder-Mead")
+ opt <- optim(fn = corphylo.LL, par = par, XX = XX, UU = UU, MM = MM, tau = tau, Vphy = Vphy, REML = REML, verbose = verbose, constrain.d = constrain.d, method = "Nelder-Mead", control = list(maxit = maxit.NM, reltol = reltol))
+
+ if (method == "SANN") {
+ opt <- optim(fn = corphylo.LL, par = par, XX = XX, UU = UU, MM = MM, tau = tau, Vphy = Vphy, REML = REML,
+ verbose = verbose, constrain.d = constrain.d, method = "SANN", control = list(maxit = maxit.SA,
+ temp = temp.SA, tmax = tmax.SA, reltol = reltol))
+ par <- opt$par
+ opt <- optim(fn = corphylo.LL, par = par, XX = XX, UU = UU, MM = MM, tau = tau, Vphy = Vphy, REML = REML,
+ verbose = verbose, constrain.d = constrain.d, method = "Nelder-Mead", control = list(maxit = maxit.NM,
+ reltol = reltol))
+ }
+
+ # Extract parameters
+ par <- Re(opt$par)
+ LL <- opt$value
+
+ L.elements <- par[1:(p + p * (p - 1)/2)]
+ L <- matrix(0, nrow = p, ncol = p)
+ L[lower.tri(L, diag = T)] <- L.elements
+ R <- t(L) %*% L
+ Rd <- diag(diag(R)^-0.5)
+ cor.matrix <- Rd %*% R %*% Rd
+
+ if (constrain.d == TRUE) {
+ logit.d <- par[(p + p * (p - 1)/2 + 1):length(par)]
+ d <- 1/(1 + exp(-logit.d))
+ } else {
+ d <- par[(p + p * (p - 1)/2 + 1):length(par)]
+ }
+
+ # OU transform
+ C <- matrix(0, nrow = p * n, ncol = p * n)
+ for (i in 1:p) for (j in 1:p) {
+ Cd <- (d[i]^tau * (d[j]^t(tau)) * (1 - (d[i] * d[j])^Vphy))/(1 - d[i] * d[j])
+ C[(n * (i - 1) + 1):(i * n), (n * (j - 1) + 1):(j * n)] <- R[i, j] * Cd
+ }
+
+ V <- C + diag(MM)
+ iV <- solve(V)
+ denom <- t(UU) %*% iV %*% UU
+ num <- t(UU) %*% iV %*% XX
+ B <- solve(denom, num)
+ B <- as.matrix(B)
+ B.cov <- solve(t(UU) %*% iV %*% UU)
+ H <- XX - UU %*% B
+
+ # Back-transform B
+ counter <- 0
+ sd.list <- matrix(0, nrow = dim(UU)[2], ncol = 1)
+ for (i in 1:p) {
+ counter <- counter + 1
+ B[counter] <- B[counter] + mean(X[, i])
+ sd.list[counter] <- sd(X[, i])
+ if (length(U) > 0) {
+ for (j in 1:ncol(U[[i]])) {
+ if (sd(U[[i]][, j]) > 0) {
+ counter <- counter + 1
+ B[counter] <- B[counter] * sd(X[, i])/sd(U[[i]][, j])
+ sd.list[counter] <- sd(X[, i])/sd(U[[i]][, j])
+ }
+ }
+ }
+ }
+ B.cov <- diag(as.numeric(sd.list)) %*% B.cov %*% diag(as.numeric(sd.list))
+ B.se <- as.matrix(diag(B.cov))^0.5
+ B.zscore <- B/B.se
+ B.pvalue <- 2 * pnorm(abs(B/B.se), lower.tail = FALSE)
+
+ # RowNames for B
+ if (length(U) > 0) {
+ B.rownames <- NULL
+ for (i in 1:p) {
+ B.rownames <- c(B.rownames, paste("B", i, ".0", sep = ""))
+ if (ncol(U[[i]]) > 0)
+ for (j in 1:ncol(U[[i]])) if (sd(U[[i]][, j]) > 0) {
+ if (is.null(colnames(U[[i]])[j]))
+ B.rownames <- c(B.rownames, paste("B", i, ".", j, sep = ""))
+ if (!is.null(colnames(U[[i]])[j]))
+ B.rownames <- c(B.rownames, paste("B", i, ".", colnames(U[[i]])[j], sep = ""))
+ }
+ }
+ } else {
+ B.rownames <- NULL
+ for (i in 1:p) {
+ B.rownames <- c(B.rownames, paste("B", i, ".0", sep = ""))
+ }
+ }
+ rownames(B) <- B.rownames
+ rownames(B.cov) <- B.rownames
+ colnames(B.cov) <- B.rownames
+ rownames(B.se) <- B.rownames
+ rownames(B.zscore) <- B.rownames
+ rownames(B.pvalue) <- B.rownames
+
+ if (REML == TRUE) {
+ logLik <- -0.5 * ((n * p) - ncol(UU)) * log(2 * pi) + 0.5 * determinant(t(XX) %*% XX)$modulus[1] - LL
+ } else {
+ logLik <- -0.5 * (n * p) * log(2 * pi) - LL
+ }
+ k <- length(par) + ncol(UU)
+ AIC <- -2 * logLik + 2 * k
+ BIC <- -2 * logLik + k * (log(n) - log(pi))
+
+ results <- list(cor.matrix = cor.matrix, d = d, B = B, B.se = B.se, B.cov = B.cov, B.zscore = B.zscore,
+ B.pvalue = B.pvalue, logLik = logLik, AIC = AIC, BIC = BIC, REML = REML, constrain.d = constrain.d,
+ XX = XX, UU = UU, MM = MM, Vphy = Vphy, R = R, V = V, C = C, convcode = opt$convergence, niter = opt$counts)
+ class(results) <- "corphylo"
+ return(results)
+}
+
+# Printing corphylo objects
+print.corphylo <- function(x, digits = max(3, getOption("digits") - 3), ...) {
+
+ cat("Call to corphylo\n\n")
+
+ logLik = x$logLik
+ AIC = x$AIC
+ BIC = x$BIC
+
+ names(logLik) = "logLik"
+ names(AIC) = "AIC"
+ names(BIC) = "BIC"
+ print(c(logLik, AIC, BIC), digits = digits)
+
+ cat("\ncorrelation matrix:\n")
+ rownames(x$cor.matrix) <- 1:dim(x$cor.matrix)[1]
+ colnames(x$cor.matrix) <- 1:dim(x$cor.matrix)[1]
+ print(x$cor.matrix, digits = digits)
+
+ cat("\nfrom OU process:\n")
+ d <- data.frame(d = x$d)
+ print(d, digits = digits)
+ if (x$constrain.d == TRUE)
+ cat("\nvalues of d constrained to be in [0, 1]\n")
+
+ cat("\ncoefficients:\n")
+ coef <- data.frame(Value = x$B, Std.Error = x$B.se, Zscore = x$B.zscore, Pvalue = x$B.pvalue)
+ rownames(coef) <- rownames(x$B)
+ printCoefmat(coef, P.values = TRUE, has.Pvalue = TRUE)
+ cat("\n")
+
+ if (x$convcode != 0)
+ cat("\nWarning: convergence in optim() not reached\n")
+}
diff --git a/R/dbd.R b/R/dbd.R
index b0beb7a..8f3c82a 100644
--- a/R/dbd.R
+++ b/R/dbd.R
@@ -1,8 +1,8 @@
-## dbd.R (2012-09-14)
+## dbd.R (2015-02-06)
## Probability Density Under Birth--Death Models
-## Copyright 2012 Emmanuel Paradis
+## Copyright 2012-2015 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
@@ -91,7 +91,10 @@ dbdTime <- function(x, birth, death, t, conditional = FALSE,
PrNt <- function(t, T, x) {
tmp <- exp(-RHO(t, T))
Wt <- tmp * (1 + INT(t))
- (1/Wt)*(1 - 1/Wt)^(x - 1)
+ out <- (1/Wt)*(1 - 1/Wt)^(x - 1)
+ zero <- x == 0
+ if (length(zero)) out[zero] <- 0
+ out
}
} else { # the unconditional case:
PrNt <- function(t, T, x)
diff --git a/R/dist.topo.R b/R/dist.topo.R
index d000da8..b50c70d 100644
--- a/R/dist.topo.R
+++ b/R/dist.topo.R
@@ -1,14 +1,46 @@
-## dist.topo.R (2014-05-16)
+## dist.topo.R (2015-04-24)
## Topological Distances, Tree Bipartitions,
## Consensus Trees, and Bootstrapping Phylogenies
-## Copyright 2005-2014 Emmanuel Paradis
+## Copyright 2005-2015 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
-dist.topo <- function(x, y, method = "PH85")
+.getTreesFromDotdotdot <- function(...)
+{
+ obj <- list(...)
+ if (length(obj) == 1 && !inherits(obj[[1]], "phylo")) obj <- obj[[1]]
+ obj
+}
+
+dist.topo <- function(x, y = NULL, method = "PH85")
+{
+ if (!is.null(y)) {
+ res <- .dist.topo(x, y, method = method)
+ } else {
+ n <- length(x) # number of trees
+ res <- numeric(n * (n - 1) /2)
+ nms <- names(x)
+ if (is.null(nms)) nms <- paste0("tree", 1:n)
+ k <- 0L
+ for (i in 1:(n - 1)) {
+ for (j in (i + 1):n) {
+ k <- k + 1L
+ res[k] <- .dist.topo(x[[i]], x[[j]], method = method)
+ }
+ }
+ attr(res, "Size") <- n
+ attr(res, "Labels") <- nms
+ attr(res, "Diag") <- attr(res, "Upper") <- FALSE
+ attr(res, "method") <- method
+ class(res) <- "dist"
+ }
+ res
+}
+
+.dist.topo <- function(x, y, method = "PH85")
{
if (method == "score" && (is.null(x$edge.length) || is.null(y$edge.length)))
stop("trees must have branch lengths for branch score distance.")
@@ -96,9 +128,7 @@ dist.topo <- function(x, y, method = "PH85")
prop.part <- function(..., check.labels = TRUE)
{
- obj <- list(...)
- if (length(obj) == 1 && !inherits(obj[[1]], "phylo")) # fix by Steve Walker (2014-03-31)
- obj <- obj[[1]]
+ obj <- .getTreesFromDotdotdot(...)
ntree <- length(obj)
if (ntree == 1) check.labels <- FALSE
if (check.labels) obj <- .compressTipLabel(obj) # fix by Klaus Schliep (2011-02-21)
@@ -133,7 +163,7 @@ print.prop.part <- function(x, ...)
summary.prop.part <- function(object, ...) attr(object, "number")
-plot.prop.part <- function(x, barcol = "blue", leftmar = 4, ...)
+plot.prop.part <- function(x, barcol = "blue", leftmar = 4, col = "red", ...)
{
if (is.null(attr(x, "labels")))
stop("cannot plot this partition object; see ?prop.part for details.")
@@ -141,23 +171,19 @@ plot.prop.part <- function(x, barcol = "blue", leftmar = 4, ...)
n <- length(attr(x, "labels"))
layout(matrix(1:2, 2, 1), heights = c(1, 3))
par(mar = c(0.1, leftmar, 0.1, 0.1))
- plot(1:L, attr(x, "number"), type = "h", col = barcol, xlim = c(1, L),
- xlab = "", ylab = "Frequency", xaxt = "n", bty = "n")
- plot(0, type = "n", xlim = c(1, L), ylim = c(1, n),
- xlab = "", ylab = "", xaxt = "n", yaxt = "n")
- for (i in 1:L) points(rep(i, length(x[[i]])), x[[i]], ...)
+ one2L <- seq_len(L)
+ plot(one2L, attr(x, "number"), type = "h", col = barcol, xlim = c(0.5, L),
+ xaxs = "i", xlab = "", ylab = "Frequency", xaxt = "n", bty = "n", ...)
+ M <- matrix(0L, L, n)
+ for (i in one2L) M[i, x[[i]]] <- 1L
+ image.default(one2L, 1:n, M, col = c("white", col), xlab = "", ylab = "", yaxt = "n")
mtext(attr(x, "labels"), side = 2, at = 1:n, las = 1)
}
prop.clades <- function(phy, ..., part = NULL, rooted = FALSE)
{
if (is.null(part)) {
- ## <FIXME>
- ## Are we going to keep the '...' way of passing trees?
- obj <- list(...)
- if (length(obj) == 1 && class(obj[[1]]) != "phylo")
- obj <- unlist(obj, recursive = FALSE)
- ## </FIXME>
+ obj <- .getTreesFromDotdotdot(...)
part <- prop.part(obj, check.labels = TRUE)
}
@@ -321,13 +347,9 @@ consensus <- function(..., p = 1, check.labels = TRUE)
pos <<- pos + size
}
}
- obj <- list(...)
- if (length(obj) == 1) {
- ## better than unlist(obj, recursive = FALSE)
- ## because "[[" keeps the class of 'obj':
- obj <- obj[[1]]
- if (class(obj) == "phylo") return(obj)
- }
+
+ obj <- .getTreesFromDotdotdot(...)
+
if (!is.null(attr(obj, "TipLabel")))
labels <- attr(obj, "TipLabel")
else {
diff --git a/R/parafit.R b/R/parafit.R
index 8c80019..a59e450 100644
--- a/R/parafit.R
+++ b/R/parafit.R
@@ -135,11 +135,11 @@ if(test.links) {
if(!is.na(stat2[k+1])) {
p.stat2 <- c(p.stat2, nGT2/(nperm2+1))
} else {
- p.stat2 <- NA
+ p.stat2 <- c(p.stat2, NA) ### Error in previous version, corrected here
}
} else {
- p.stat1 <- NA
- p.stat2 <- NA
+ p.stat1 <- c(p.stat1, NA) ### Error in previous version, corrected here
+ p.stat2 <- c(p.stat2, NA) ### Error in previous version, corrected here
}
}
#
diff --git a/R/plot.phylo.R b/R/plot.phylo.R
index 9c17ca0..fd6f75e 100644
--- a/R/plot.phylo.R
+++ b/R/plot.phylo.R
@@ -1,8 +1,8 @@
-## plot.phylo.R (2014-08-21)
+## plot.phylo.R (2015-02-25)
## Plot Phylogenies
-## Copyright 2002-2014 Emmanuel Paradis
+## Copyright 2002-2015 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
@@ -56,10 +56,6 @@ plot.phylo <-
## the order of the last two conditions is important:
if (type %in% c("unrooted", "radial") || !use.edge.length ||
is.null(x$root.edge) || !x$root.edge) root.edge <- FALSE
- if (type == "fan" && root.edge) {
- warning("drawing root edge with type = 'fan' is not yet supported")
- root.edge <- FALSE
- }
phyloORclado <- type %in% c("phylogram", "cladogram")
horizontal <- direction %in% c("rightwards", "leftwards")
@@ -145,6 +141,7 @@ if (phyloORclado) {
r <- 1/r
}
theta <- theta + rotate.tree
+ if (root.edge) r <- r + x$root.edge
xx <- r * cos(theta)
yy <- r * sin(theta)
}, "unrooted" = {
@@ -334,12 +331,19 @@ if (plot) {
} else
cladogram.plot(x$edge, xx, yy, edge.color, edge.width, edge.lty)
}
- if (root.edge)
- switch(direction,
- "rightwards" = segments(0, yy[ROOT], x$root.edge, yy[ROOT]),
- "leftwards" = segments(xx[ROOT], yy[ROOT], xx[ROOT] + x$root.edge, yy[ROOT]),
- "upwards" = segments(xx[ROOT], 0, xx[ROOT], x$root.edge),
- "downwards" = segments(xx[ROOT], yy[ROOT], xx[ROOT], yy[ROOT] + x$root.edge))
+ if (root.edge) {
+ rootcol <- if (length(edge.color) == 1) edge.color else "black"
+ if (type == "fan") {
+ tmp <- polar2rect(x$root.edge, theta[ROOT])
+ segments(0, 0, tmp$x, tmp$y, col = rootcol)
+ } else {
+ switch(direction,
+ "rightwards" = segments(0, yy[ROOT], x$root.edge, yy[ROOT], col = rootcol),
+ "leftwards" = segments(xx[ROOT], yy[ROOT], xx[ROOT] + x$root.edge, yy[ROOT], col = rootcol),
+ "upwards" = segments(xx[ROOT], 0, xx[ROOT], x$root.edge, col = rootcol),
+ "downwards" = segments(xx[ROOT], yy[ROOT], xx[ROOT], yy[ROOT] + x$root.edge, col = rootcol))
+ }
+ }
if (show.tip.label) {
if (is.expression(x$tip.label)) underscore <- TRUE
if (!underscore) x$tip.label <- gsub("_", " ", x$tip.label)
@@ -670,20 +674,14 @@ node.height <- function(phy, clado.style = FALSE)
as.double(yy))[[6]]
}
-node.height.clado <- function(phy)
-{
- warning("the function 'node.height.clado' will be removed soon.\nUse node.height(phy, clado.style = TRUE) instead.")
- node.height(phy, TRUE)
-}
-
plot.multiPhylo <- function(x, layout = 1, ...)
{
layout(matrix(1:layout, ceiling(sqrt(layout)), byrow = TRUE))
- if (!devAskNewPage() && interactive()) {
+ if (!devAskNewPage() && !names(dev.cur()) %in% c("pdf", "postscript")) {
devAskNewPage(TRUE)
on.exit(devAskNewPage(FALSE))
}
- for (i in 1:length(x)) plot(x[[i]], ...)
+ for (i in seq_along(x)) plot(x[[i]], ...)
}
trex <- function(phy, title = TRUE, subbg = "lightyellow3",
diff --git a/R/read.tree.R b/R/read.tree.R
index 391ed73..4ec8454 100644
--- a/R/read.tree.R
+++ b/R/read.tree.R
@@ -1,4 +1,4 @@
-## read.tree.R (2012-09-14)
+## read.tree.R (2015-01-12)
## Read Tree Files in Parenthetic Format
@@ -133,8 +133,10 @@ read.tree <- function(file = "", text = NULL, tree.names = NULL, skip = 0,
## Suggestion from Olivier Francois (added 2006-07-15):
if (is.na(y[1])) return(NULL)
STRING <- character(Ntree)
- for (i in 1:Ntree)
- STRING[i] <- paste(tree[x[i]:y[i]], sep = "", collapse = "")
+ for (i in 1:Ntree) {
+ tmp <- paste(tree[x[i]:y[i]], sep = "", collapse = "")
+ STRING[i] <- gsub("\\[[^]]*\\]", "", tmp) # delete comments (fix 2015-01-12)
+ }
tmp <- unlist(lapply(STRING, unname))
tmpnames <- tmp[c(TRUE, FALSE)]
diff --git a/R/root.R b/R/root.R
index 681d1de..d0eef6f 100644
--- a/R/root.R
+++ b/R/root.R
@@ -1,8 +1,8 @@
-## root.R (2014-05-28)
+## root.R (2015-04-08)
## Root of Phylogenetic Trees
-## Copyright 2004-2014 Emmanuel Paradis
+## Copyright 2004-2015 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
@@ -144,6 +144,10 @@ root <- function(phy, outgroup, node = NULL,
if (newroot == ROOT) {
## assumes length(outgroup) == 1
if (resolve.root) {
+ ## add this check 2015-04-08:
+ if (!is.null(node))
+ stop("ambiguous resolution of the root node: please specify an explicit outgroup")
+ ##
snw <- which(phy$edge[, 1L] == newroot)
if (length(snw) > 2) {
i <- which(phy$edge[, 2L] == outgroup) # see comment above
diff --git a/R/write.nexus.R b/R/write.nexus.R
index ac383f7..0d6b61c 100644
--- a/R/write.nexus.R
+++ b/R/write.nexus.R
@@ -1,23 +1,16 @@
-## write.nexus.R (2012-03-30)
+## write.nexus.R (2014-12-10)
## Write Tree File in Nexus Format
-## Copyright 2003-2012 Emmanuel Paradis
+## Copyright 2003-2014 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
write.nexus <- function(..., file = "", translate = TRUE)
{
- obj <- list(...)
- ## We insure that all trees are in a list, even if there is a single one:
- if (length(obj) == 1) {
- if (class(obj[[1]]) == "phylo") ntree <- 1
- else {
- obj <- obj[[1]] # NOT use unlist()
- ntree <- length(obj)
- }
- } else ntree <- length(obj)
+ obj <- .getTreesFromDotdotdot(...)
+ ntree <- length(obj)
cat("#NEXUS\n", file = file)
cat(paste("[R-package APE, ", date(), "]\n\n", sep = ""),
file = file, append = TRUE)
diff --git a/R/write.nexus.data.R b/R/write.nexus.data.R
index 0223c76..2be1dd4 100644
--- a/R/write.nexus.data.R
+++ b/R/write.nexus.data.R
@@ -1,8 +1,8 @@
-## write.nexus.data.R (2014-09-24)
+## write.nexus.data.R (2015-02-04)
## Write Character Data in NEXUS Format
-## Copyright 2006-2014 Johan Nylander, Emmanuel Paradis
+## Copyright 2006-2015 Johan Nylander, Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
@@ -77,7 +77,7 @@ write.nexus.data <-
NCHAR <- paste("NCHAR=", nchars, sep = "")
NTAX <- paste0("NTAX=", ntax)
- DATATYPE <- paste0("DATATYPE", format)
+ DATATYPE <- paste0("DATATYPE=", format) # fix by Robin Cristofari (2015-02-04)
if (is.null(charsperline)) {
if (nchars <= defcharsperline) {
diff --git a/build/vignette.rds b/build/vignette.rds
index a13ffe5..dd76fc1 100644
Binary files a/build/vignette.rds and b/build/vignette.rds differ
diff --git a/data/hivtree.newick.rda b/data/hivtree.newick.rda
index c68fd4b..1662abc 100644
Binary files a/data/hivtree.newick.rda and b/data/hivtree.newick.rda differ
diff --git a/data/landplants.newick.rda b/data/landplants.newick.rda
index 9ef4514..1f40b47 100644
Binary files a/data/landplants.newick.rda and b/data/landplants.newick.rda differ
diff --git a/data/opsin.newick.rda b/data/opsin.newick.rda
index eec50af..e7089dc 100644
Binary files a/data/opsin.newick.rda and b/data/opsin.newick.rda differ
diff --git a/inst/doc/MoranI.pdf b/inst/doc/MoranI.pdf
index 15a3d25..4e5b25e 100644
Binary files a/inst/doc/MoranI.pdf and b/inst/doc/MoranI.pdf differ
diff --git a/man/DNAbin.Rd b/man/DNAbin.Rd
index 020f6f7..1cbce70 100644
--- a/man/DNAbin.Rd
+++ b/man/DNAbin.Rd
@@ -80,7 +80,7 @@
}
\references{
Paradis, E. (2007) A Bit-Level Coding Scheme for Nucleotides.
- \url{http://ape.mpl.ird.fr/misc/BitLevelCodingScheme_20April2007.pdf}
+ \url{http://ape-package.ird.fr/misc/BitLevelCodingScheme_20April2007.pdf}
Paradis, E. (2012) \emph{Analysis of Phylogenetics and Evolution with
R (Second Edition).} New York: Springer.
diff --git a/man/LTT.Rd b/man/LTT.Rd
index 7ef9012..7d76b17 100644
--- a/man/LTT.Rd
+++ b/man/LTT.Rd
@@ -38,9 +38,13 @@ LTT(birth = 0.1, death = 0, N = 100, Tmax = 50, PI = 95,
\details{
For the moment, this works well when \code{birth} and \code{death} are
constant. Some improvements are under progress for time-dependent
- rates (but see below for and example).
+ rates (but see below for an example).
}
\references{
+ Hallinan, N. (2012) The generalized time variable reconstructed
+ birth--death process. \emph{Journal of Theoretical Biology},
+ \bold{300}, 265--276.
+
Paradis, E. (2011) Time-dependent speciation and extinction from
phylogenies: a least squares approach. \emph{Evolution}, \bold{65},
661--672.
diff --git a/man/MoranI.Rd b/man/MoranI.Rd
index fd7b28e..adcc998 100644
--- a/man/MoranI.Rd
+++ b/man/MoranI.Rd
@@ -60,9 +60,7 @@
}
\author{Julien Dutheil \email{julien.dutheil at univ-montp2.fr} and
Emmanuel Paradis}
-\seealso{
- \code{\link{weight.taxo}}
-}
+\seealso{\code{\link{weight.taxo}}}
\examples{
tr <- rtree(30)
x <- rnorm(30)
diff --git a/man/ace.Rd b/man/ace.Rd
index c61b7f6..7094091 100644
--- a/man/ace.Rd
+++ b/man/ace.Rd
@@ -127,8 +127,8 @@ ace(x, phy, type = "continuous", method = if (type == "continuous")
discrete characters are computed with a joint estimation procedure
using a procedure similar to the one described in Pupko et al. (2000).
If \code{marginal = TRUE}, a marginal estimation procedure is used
- (this was the only choice until ape 3.1-1). With this second method,
- the likelihood values at a given node are computed using only the
+ (this was the only choice until ape 3.1-1). With this method, the
+ likelihood values at a given node are computed using only the
information from the tips (and branches) descending from this node.
With the joint estimation, all information is used for each node. The
difference between these two methods is further explained in
@@ -138,7 +138,7 @@ ace(x, phy, type = "continuous", method = if (type == "continuous")
estimates of both methods are very close.
With discrete characters it is necessary to compute the exponential of
- the rate matrix. The only possibility until \pkg{ape} 3.0-7) was the
+ the rate matrix. The only possibility until \pkg{ape} 3.0-7 was the
function \code{\link{matexpo}} in \pkg{ape}. If \code{use.expm = TRUE}
and \code{use.eigen = FALSE}, the function \code{\link[expm]{expm}},
in the package of the same name, is used. \code{matexpo} is faster but
@@ -148,6 +148,14 @@ ace(x, phy, type = "continuous", method = if (type == "continuous")
exponential; see details in Lebl (2013, sect. 3.8.3). This is much
faster and is now the default.
}
+\note{
+ Liam Revell points out that for discrete characters the ancestral
+ likelihood values returned with \code{marginal = FALSE} are actually
+ the marginal estimates, while setting \code{marginal = TRUE} returns
+ the conditional (scaled) likelihoods of the subtree:
+
+ \url{http://blog.phytools.org/2015/05/about-how-acemarginaltrue-does-not.html}
+}
\value{
an object of class \code{"ace"} with the following elements:
diff --git a/man/as.bitsplits.Rd b/man/as.bitsplits.Rd
index 4aeadb8..d9cb25e 100644
--- a/man/as.bitsplits.Rd
+++ b/man/as.bitsplits.Rd
@@ -2,6 +2,7 @@
\alias{as.bitsplits}
\alias{as.bitsplits.prop.part}
\alias{print.bitsplits}
+\alias{sort.bitsplits}
\alias{bitsplits}
\alias{countBipartitions}
\alias{as.prop.part}
@@ -24,6 +25,7 @@ countBipartitions(phy, X)
as.bitsplits(x)
\method{as.bitsplits}{prop.part}(x)
\method{print}{bitsplits}(x, ...)
+\method{sort}{bitsplits}(x, decreasing = FALSE, ...)
as.prop.part(x)
\method{as.prop.part}{bitsplits}(x)
}
@@ -31,6 +33,8 @@ as.prop.part(x)
\item{x}{an object of the appropriate class.}
\item{phy}{an object of class \code{"phylo"}.}
\item{X}{an object of class \code{"multiPhylo"}.}
+ \item{decreasing}{a logical value to sort the bipartitions in
+ increasing (the default) or decreasing order of their frequency.}
\item{\dots}{further arguments passed to or from other methods.}
}
\details{
@@ -40,8 +44,8 @@ as.prop.part(x)
on ape's web site.
}
\value{
- \code{bitsplits} and \code{as.bitsplits} return an object of class
- \code{"bitsplits"}.
+ \code{bitsplits}, \code{as.bitsplits}, and \code{sort} return an object
+ of class \code{"bitsplits"}.
\code{countBipartitions} returns a vector of integers.
diff --git a/man/binaryPGLMM.Rd b/man/binaryPGLMM.Rd
new file mode 100644
index 0000000..6a151a9
--- /dev/null
+++ b/man/binaryPGLMM.Rd
@@ -0,0 +1,268 @@
+\name{binaryPGLMM}
+\alias{binaryPGLMM}
+\alias{binaryPGLMM.sim}
+\alias{print.binaryPGLMM}
+\title{Phylogenetic Generalized Linear Mixed Model for Binary Data}
+\description{
+ binaryPGLMM performs linear regression for binary phylogenetic data, estimating regression coefficients with approximate standard errors. It simultaneously estimates the strength of phylogenetic signal in the residuals and gives an approximate conditional likelihood ratio test for the hypothesis that there is no signal. Therefore, when applied without predictor (independent) variables, it gives a test for phylogenetic signal for binary data. The method uses a GLMM approach, alternating [...]
+
+ binaryPGLMM.sim is a companion function that simulates binary phylogenetic data of the same structure analyzed by binaryPGLMM.
+}
+\usage{
+binaryPGLMM(formula, data = list(), phy, s2.init = 0.1,
+ B.init = NULL, tol.pql = 10^-6, maxit.pql = 200,
+ maxit.reml = 100)
+
+binaryPGLMM.sim(formula, data = list(), phy, s2 = NULL, B = NULL, nrep = 1)
+
+\method{print}{binaryPGLMM}(x, digits = max(3, getOption("digits") - 3), ...)
+}
+\arguments{
+ \item{formula}{a two-sided linear formula object describing the
+ fixed-effects of the model; for example, Y ~ X.}
+ \item{data}{a data frame containing the variables named in formula.}
+ \item{phy}{a phylogenetic tree as an object of class "phylo".}
+ \item{s2.init}{an initial estimate of s2, the scaling component of the
+ variance in the PGLMM. A value of s2 = 0 implies no phylogenetic
+ signal. Note that the variance-covariance matrix given by the
+ phylogeny phy is scaled to have determinant = 1.}
+ \item{B.init}{initial estimates of B, the matrix containing regression
+ coefficients in the model. This matrix must have
+ dim(B.init)=c(p+1,1), where p is the number of predictor
+ (independent) variables; the first element of B corresponds to the
+ intercept, and the remaining elements correspond in order to the
+ predictor (independent) variables in the model.}
+ \item{tol.pql}{a control parameter dictating the tolerance for
+ convergence for the PQL optimization.}
+ \item{maxit.pql}{a control parameter dictating the maximum number of
+ iterations for the PQL optimization.}
+ \item{maxit.reml}{a control parameter dictating the maximum number of
+ iterations for the REML optimization.}
+ \item{x}{an object of class "binaryPGLMM".}
+ \item{s2}{in binaryPGLMM.sim, value of s2. See s2.init.}
+ \item{B}{in binaryPGLMM.sim, value of B, the matrix containing regression coefficients in the model. See B.init.}
+ \item{nrep}{in binaryPGLMM.sim, number of compete data sets produced.}
+ \item{digits}{the number of digits to print.}
+ \item{\dots}{further arguments passed to \code{print}.}
+}
+\details{
+ The function estimates parameters for the model
+
+ \deqn{Pr(Y = 1) = q }
+ \deqn{q = inverse.logit(b0 + b1 * x1 + b2 * x2 + \dots + \epsilon)}
+ \deqn{\epsilon ~ Gaussian(0, s2 * V) }
+
+where \eqn{V} is a variance-covariance matrix derived from a phylogeny (typically under the assumption of Brownian motion evolution). Although mathematically there is no requirement for \eqn{V} to be ultrametric, forcing \eqn{V} into ultrametric form can aide in the interpretation of the model, because in regression for binary dependent variables, only the off-diagonal elements (i.e., covariances) of matrix \eqn{V} are biologically meaningful (see Ives & Garland 2014).
+
+The function converts a phylo tree object into a variance-covariance matrix, and further standardizes this matrix to have determinant = 1. This in effect standardizes the interpretation of the scalar s2. Although mathematically not required, it is a very good idea to standardize the predictor (independent) variables to have mean 0 and variance 1. This will make the function more robust and improve the interpretation of the regression coefficients. For categorical (factor) predictor varia [...]
+
+The estimation method alternates between PQL to obtain estimates of the mean components of the model (this is the standard approach to estimating GLMs) and REML to obtain estimates of the variance components. This method gives relatively fast and robust estimation. Nonetheless, the estimates of the coefficients B will generally be upwards bias, as is typical of estimation for binary data. The standard errors of B are computed from the PQL results conditional on the estimate of s2 and the [...]
+
+It is a good idea to confirm statistical inferences using parametric bootstrapping, and the companion function binaryPGLMM.sim gives a simply tool for this. See Examples below.
+}
+\value{
+ An object of class "binaryPGLMM".
+
+ \item{formula}{formula specifying the regression model.}
+ \item{B}{estimates of the regression coefficients.}
+ \item{B.se}{approximate PQL standard errors of the regression
+ coefficients.}
+ \item{B.cov}{approximate PQL covariance matrix for the regression
+ coefficients.}
+ \item{B.zscore}{approximate PQL Z scores for the regression
+ coefficients.}
+ \item{B.pvalue}{approximate PQL tests for the regression coefficients
+ being different from zero.}
+ \item{s2}{phylogenetic signal measured as the scalar magnitude of the
+ phylogenetic variance-covariance matrix s2 * V.}
+ \item{P.H0.s2}{approximate likelihood ratio test of the hypothesis H0
+ that s2 = 0. This test is based on the conditional REML (keeping the
+ regression coefficients fixed) and is prone to inflated type 1 errors.}
+ \item{mu}{for each data point y, the estimate of p that y = 1.}
+ \item{b}{for each data point y, the estimate of inverse.logit(p).}
+ \item{X}{the predictor (independent) variables returned in matrix form
+ (including 1s in the first column).}
+ \item{H}{residuals of the form b + (Y - mu)/(mu * (1 - mu)).}
+ \item{B.init}{the user-provided initial estimates of B. If B.init is
+ not provided, these are estimated using glm() assuming no phylogenetic
+ signal. The glm() estimates can generate convergence problems, so
+ using small values (e.g., 0.01) is more robust but slower.}
+ \item{VCV}{the standardized phylogenetic variance-covariance matrix.}
+ \item{V}{estimate of the covariance matrix of H.}
+ \item{convergeflag}{flag for cases when convergence failed.}
+ \item{iteration}{number of total iterations performed.}
+ \item{converge.test.B}{final tolerance for B.}
+ \item{converge.test.s2}{final tolerance for s2.}
+ \item{rcondflag}{number of times B is reset to 0.01. This is done when
+ rcond(V) < 10^(-10), which implies that V cannot be inverted.}
+ \item{Y}{in binaryPGLMM.sim, the simulated values of Y.}
+}
+\author{Anthony R. Ives}
+\references{
+ Ives, A. R. and Helmus, M. R. (2011) Generalized linear mixed models
+ for phylogenetic analyses of community structure. \emph{Ecological
+ Monographs}, \bold{81}, 511--525.
+
+ Ives, A. R. and Garland, T., Jr. (2014) Phylogenetic regression for
+ binary dependent variables. Pages 231--261 \emph{in} L. Z. Garamszegi,
+ editor. \emph{Modern Phylogenetic Comparative Methods and Their
+ Application in Evolutionary Biology}. Springer-Verlag, Berlin
+ Heidelberg.
+
+}
+\seealso{
+ package \pkg{pez} and its function \code{communityPGLMM};
+ package \pkg{phylolm} and its function \code{phyloglm};
+ package \pkg{MCMCglmm}
+}
+\examples{
+## Illustration of binaryPGLMM() with simulated data
+
+# Generate random phylogeny
+
+n <- 100
+phy <- compute.brlen(rtree(n=n), method = "Grafen", power = 1)
+
+# Generate random data and standardize to have mean 0 and variance 1
+X1 <- rTraitCont(phy, model = "BM", sigma = 1)
+X1 <- (X1 - mean(X1))/var(X1)
+
+# Simulate binary Y
+sim.dat <- data.frame(Y=array(0, dim=n), X1=X1, row.names=phy$tip.label)
+sim.dat$Y <- binaryPGLMM.sim(Y ~ X1, phy=phy, data=sim.dat, s2=.5,
+ B=matrix(c(0,.25),nrow=2,ncol=1), nrep=1)$Y
+
+# Fit model
+binaryPGLMM(Y ~ X1, phy=phy, data=sim.dat)
+
+\dontrun{
+# Compare with phyloglm()
+library(phylolm)
+summary(phyloglm(Y ~ X1, phy=phy, data=sim.dat))
+
+# Compare with glm() that does not account for phylogeny
+summary(glm(Y ~ X1, data=sim.dat, family="binomial"))
+
+# Compare with logistf() that does not account
+# for phylogeny but is less biased than glm()
+library(logistf)
+logistf(Y ~ X1, data=sim.dat)
+
+# Compare with MCMCglmm
+library(MCMCglmm)
+
+V <- vcv(phy)
+V <- V/max(V)
+detV <- exp(determinant(V)$modulus[1])
+V <- V/detV^(1/n)
+
+invV <- Matrix(solve(V),sparse=T)
+sim.dat$species <- phy$tip.label
+rownames(invV) <- sim.dat$species
+
+nitt <- 43000
+thin <- 10
+burnin <- 3000
+
+prior <- list(R=list(V=1, fix=1), G=list(G1=list(V=1, nu=1000, alpha.mu=0, alpha.V=1)))
+summary(MCMCglmm(Y ~ X1, random=~species, ginvers=list(species=invV),
+ data=sim.dat, slice=TRUE, nitt=nitt, thin=thin, burnin=burnin,
+ family="categorical", prior=prior, verbose=FALSE))
+
+## Examine bias in estimates of B1 and s2 from binaryPGLMM with
+# simulated data. Note that this will take a while.
+
+Reps = 1000
+
+s2 <- 0.4
+B1 <- 1
+
+meanEsts <- data.frame(n = Inf, B1 = B1, s2 = s2, Pr.s2 = 1, propconverged = 1)
+
+for (n in c(160, 80, 40, 20)) {
+
+ meanEsts.n <- data.frame(B1 = 0, s2 = 0, Pr.s2 = 0, convergefailure = 0)
+ for (rep in 1:Reps) {
+ phy <- compute.brlen(rtree(n = n), method = "Grafen", power = 1)
+ X <- rTraitCont(phy, model = "BM", sigma = 1)
+ X <- (X - mean(X))/var(X)
+
+ sim.dat <- data.frame(Y = array(0, dim = n), X = X, row.names = phy$tip.label)
+ sim <- binaryPGLMM.sim(Y ~ 1 + X, phy = phy, data = sim.dat, s2 = s2,
+ B = matrix(c(0,B1), nrow = 2, ncol = 1), nrep = 1)
+ sim.dat$Y <- sim$Y
+
+ z <- binaryPGLMM(Y ~ 1 + X, phy = phy, data = sim.dat)
+
+ meanEsts.n[rep, ] <- c(z$B[2], z$s2, z$P.H0.s2, z$convergeflag == "converged")
+ }
+converged <- meanEsts.n[,4]
+meanEsts <- rbind(meanEsts,
+ c(n, mean(meanEsts.n[converged==1,1]),
+ mean(meanEsts.n[converged==1,2]),
+ mean(meanEsts.n[converged==1, 3] < 0.05),
+ mean(converged)))
+}
+meanEsts
+
+# Results output for B1 = 0.5, s2 = 0.4; n-Inf gives the values used to
+# simulate the data
+# n B1 s2 Pr.s2 propconverged
+# 1 Inf 1.000000 0.4000000 1.00000000 1.000
+# 2 160 1.012719 0.4479946 0.36153072 0.993
+# 3 80 1.030876 0.5992027 0.24623116 0.995
+# 4 40 1.110201 0.7425203 0.13373860 0.987
+# 5 20 1.249886 0.8774708 0.05727377 0.873
+
+
+## Examine type I errors for estimates of B0 and s2 from binaryPGLMM()
+# with simulated data. Note that this will take a while.
+
+Reps = 1000
+
+s2 <- 0
+B0 <- 0
+B1 <- 0
+
+H0.tests <- data.frame(n = Inf, B0 = B0, s2 = s2, Pr.B0 = .05,
+ Pr.s2 = .05, propconverged = 1)
+for (n in c(160, 80, 40, 20)) {
+
+ ests.n <- data.frame(B1 = 0, s2 = 0, Pr.B0 = 0, Pr.s2 = 0, convergefailure = 0)
+ for (rep in 1:Reps) {
+ phy <- compute.brlen(rtree(n = n), method = "Grafen", power = 1)
+ X <- rTraitCont(phy, model = "BM", sigma = 1)
+ X <- (X - mean(X))/var(X)
+
+ sim.dat <- data.frame(Y = array(0, dim = n), X = X, row.names = phy$tip.label)
+ sim <- binaryPGLMM.sim(Y ~ 1, phy = phy, data = sim.dat, s2 = s2,
+ B = matrix(B0, nrow = 1, ncol = 1), nrep = 1)
+ sim.dat$Y <- sim$Y
+
+ z <- binaryPGLMM(Y ~ 1, phy = phy, data = sim.dat)
+
+ ests.n[rep, ] <- c(z$B[1], z$s2, z$B.pvalue, z$P.H0.s2, z$convergeflag == "converged")
+ }
+
+converged <- ests.n[,5]
+H0.tests <- rbind(H0.tests,
+ c(n, mean(ests.n[converged==1,1]),
+ mean(ests.n[converged==1,2]),
+ mean(ests.n[converged==1, 3] < 0.05),
+ mean(ests.n[converged==1, 4] < 0.05),
+ mean(converged)))
+}
+H0.tests
+
+# Results for type I errors for B0 = 0 and s2 = 0; n-Inf gives the values
+# used to simulate the data. These results show that binaryPGLMM() tends to
+# have lower-than-nominal p-values; fewer than 0.05 of the simulated
+# data sets have H0:B0=0 and H0:s2=0 rejected at the alpha=0.05 level.
+# n B0 s2 Pr.B0 Pr.s2 propconverged
+# 1 Inf 0.0000000000 0.00000000 0.05000000 0.05000000 1.000
+# 2 160 -0.0009350357 0.07273163 0.02802803 0.04804805 0.999
+# 3 80 -0.0085831477 0.12205876 0.04004004 0.03403403 0.999
+# 4 40 0.0019303847 0.25486307 0.02206620 0.03711133 0.997
+# 5 20 0.0181394905 0.45949266 0.02811245 0.03313253 0.996
+}}
+\keyword{regression}
diff --git a/man/boot.phylo.Rd b/man/boot.phylo.Rd
index b031e0a..8c699f4 100644
--- a/man/boot.phylo.Rd
+++ b/man/boot.phylo.Rd
@@ -27,7 +27,7 @@ prop.part(..., check.labels = TRUE)
prop.clades(phy, ..., part = NULL, rooted = FALSE)
\method{print}{prop.part}(x, ...)
\method{summary}{prop.part}(object, ...)
-\method{plot}{prop.part}(x, barcol = "blue", leftmar = 4, ...)
+\method{plot}{prop.part}(x, barcol = "blue", leftmar = 4, col = "red", ...)
}
\arguments{
\item{phy}{an object of class \code{"phylo"}.}
@@ -57,6 +57,7 @@ prop.clades(phy, ..., part = NULL, rooted = FALSE)
partitions in the upper panel.}
\item{leftmar}{the size of the margin on the left to display the tip
labels.}
+ \item{col}{the colour used to visualise the bipartitions.}
}
\details{
The argument \code{FUN} in \code{boot.phylo} must be the function used
diff --git a/man/coalescent.intervals.Rd b/man/coalescent.intervals.Rd
index 1a3ac19..9ea0b14 100644
--- a/man/coalescent.intervals.Rd
+++ b/man/coalescent.intervals.Rd
@@ -2,7 +2,6 @@
\alias{coalescent.intervals}
\alias{coalescent.intervals.phylo}
\alias{coalescent.intervals.default}
-
\title{Coalescent Intervals}
\usage{
coalescent.intervals(x)
@@ -33,16 +32,12 @@ An object of class \code{"coalescentIntervals"} with the following entries:
\code{\link{branching.times}}, \code{\link{collapsed.intervals}},
\code{\link{read.tree}}.
}
-
-\author{Korbinian Strimmer (\url{http://www.stat.uni-muenchen.de/~strimmer/})}
-
+\author{Korbinian Strimmer}
\examples{
data("hivtree.newick") # example tree in NH format
tree.hiv <- read.tree(text = hivtree.newick) # load tree
-
ci <- coalescent.intervals(tree.hiv) # from tree
ci
-
data("hivtree.table") # same tree, but in table format
ci <- coalescent.intervals(hivtree.table$size) # from vector of interval lengths
ci
diff --git a/man/collapsed.intervals.Rd b/man/collapsed.intervals.Rd
index 27a79ce..1b0e1d6 100644
--- a/man/collapsed.intervals.Rd
+++ b/man/collapsed.intervals.Rd
@@ -1,6 +1,5 @@
\name{collapsed.intervals}
\alias{collapsed.intervals}
-
\title{Collapsed Coalescent Intervals}
\usage{
collapsed.intervals(ci, epsilon=0)
@@ -25,7 +24,6 @@ interval is larger than \code{epsilon}. Note that this approach prevents the
occurrence of zero-length intervals at the present.
For more details see Strimmer and Pybus (2001).
}
-
\value{
An object of class \code{"collapsedIntervals"} with the following entries:
@@ -42,33 +40,24 @@ An object of class \code{"collapsedIntervals"} with the following entries:
intervals.}
\item{epsilon}{The value of the underlying smoothing parameter.}
}
-
-\author{Korbinian Strimmer (\url{http://www.stat.uni-muenchen.de/~strimmer/})}
-
+\author{Korbinian Strimmer}
\seealso{
\code{\link{coalescent.intervals}},\code{\link{skyline}}.
}
-
-
\references{
Strimmer, K. and Pybus, O. G. (2001) Exploring the demographic history
of DNA sequences using the generalized skyline plot. \emph{Molecular
Biology and Evolution}, \bold{18}, 2298--2305.
}
-
\examples{
data("hivtree.table") # example tree
-
# colescent intervals from vector of interval lengths
ci <- coalescent.intervals(hivtree.table$size)
ci
-
# collapsed intervals
cl1 <- collapsed.intervals(ci,0)
cl2 <- collapsed.intervals(ci,0.0119)
-
cl1
cl2
-
}
\keyword{manip}
diff --git a/man/corphylo.Rd b/man/corphylo.Rd
new file mode 100644
index 0000000..5ea27ff
--- /dev/null
+++ b/man/corphylo.Rd
@@ -0,0 +1,224 @@
+\name{corphylo}
+\alias{corphylo}
+\alias{print.corphylo}
+\title{Correlations among Multiple Traits with Phylogenetic Signal}
+\description{
+ This function calculates Pearson correlation coefficients for multiple continuous traits that may have phylogenetic signal, allowing users to specify measurement error as the standard error of trait values at the tips of the phylogenetic tree. Phylogenetic signal for each trait is estimated from the data assuming that trait evolution is given by a Ornstein-Uhlenbeck process. Thus, the function allows the estimation of phylogenetic signal in multiple traits while incorporating correla [...]
+}
+\usage{
+corphylo(X, U = list(), SeM = NULL, phy = NULL, REML = TRUE,
+method = c("Nelder-Mead", "SANN"), constrain.d = FALSE, reltol = 10^-6,
+maxit.NM = 1000, maxit.SA = 1000, temp.SA = 1, tmax.SA = 1, verbose = FALSE)
+
+\method{print}{corphylo}(x, digits = max(3, getOption("digits") - 3), ...)
+}
+\arguments{
+ \item{X}{a n x p matrix with p columns containing the values for the n taxa. Rows of X should have rownames matching the taxon names in phy.}
+ \item{U}{a list of p matrices corresponding to the p columns of X, with each matrix containing independent variables for the corresponding column of X. The rownames of each matrix within U must be the same as X, or alternatively, the order of values in rows must match those in X. If U is omitted, only the mean (aka intercept) for each column of X is estimated. If U[[i]] is NULL, only an intercept is estimated for X[, i]. If all values of U[[i]][j] are the same, this variable is automat [...]
+ \item{SeM}{a n x p matrix with p columns containing standard errors of the trait values in X. The rownames of SeM must be the same as X, or alternatively, the order of values in rows must match those in X. If SeM is omitted, the trait values are assumed to be known without error. If only some traits have mesurement errors, the remaining traits can be given zero-valued standard errors.}
+ \item{phy}{a phylo object giving the phylogenetic tree. The rownames of phy must be the same as X, or alternatively, the order of values in rows must match those in X.}
+ \item{REML}{whether REML or ML is used for model fitting.}
+ \item{method}{in optim(), either Nelder-Mead simplex minimization or SANN (simulated annealing) minimization is used. If SANN is used, it is followed by Nelder-Mead minimization.}
+ \item{constrain.d}{if constrain.d is TRUE, the estimates of d are constrained to be between zero and 1. This can make estimation more stable and can be tried if convergence is problematic. This does not necessarily lead to loss of generality of the results, because before using corphylo, branch lengths of phy can be transformed so that the "starter" tree has strong phylogenetic signal.}
+ \item{reltol}{a control parameter dictating the relative tolerance for convergence in the optimization; see optim().}
+ \item{maxit.NM}{a control parameter dictating the maximum number of iterations in the optimization with Nelder-Mead minimization; see optim().}
+ \item{maxit.SA}{a control parameter dictating the maximum number of iterations in the optimization with SANN minimization; see optim().}
+ \item{temp.SA}{a control parameter dictating the starting temperature in the optimization with SANN minimization; see optim().}
+ \item{tmax.SA}{a control parameter dictating the number of function evaluations at each temperature in the optimization with SANN minimization; see optim().}
+ \item{verbose}{if TRUE, the model logLik and running estimates of the
+ correlation coefficients and values of d are printed each iteration
+ during optimization.}
+ \item{x}{an objects of class corphylo.}
+ \item{digits}{the number of digits to be printed.}
+ \item{\dots}{arguments passed to and from other methods.}
+}
+\details{
+ For the case of two variables, the function estimates parameters for the model of the form, for example,
+
+ \deqn{X[1] = B[1,0] + B[1,1] * u[1,1] + \epsilon[1]}
+ \deqn{X[2] = B[2,0] + B[2,1] * u[2,1] + \epsilon[2]}
+ \deqn{\epsilon ~ Gaussian(0, V) }
+
+where \eqn{B[1,0]}, \eqn{B[1,1]}, \eqn{B[2,0]}, and \eqn{B[2,1]} are regression coefficients, and \eqn{V} is a variance-covariance matrix containing the correlation coefficient r, parameters of the OU process \eqn{d1} and \eqn{d2}, and diagonal matrices \eqn{M1} and \eqn{M2} of measurement standard errors for \eqn{X[1]} and \eqn{X[2]}. The matrix \eqn{V} is \eqn{2n x 2n}, with \eqn{n x n} blocks given by
+
+\deqn{V[1,1] = C[1,1](d1) + M1}
+\deqn{V[1,2] = C[1,2](d1,d2)}
+\deqn{V[2,1] = C[2,1](d1,d2)}
+\deqn{V[2,2] = C[2,2](d2) + M2}
+
+where \eqn{C[i,j](d1,d2)} are derived from phy under the assumption of joint OU evolutionary processes for each trait (see Zheng et al. 2009). This formulation extends in the obvious way to more than two traits.
+
+}
+\value{
+ An object of class "corphylo".
+
+ \item{cor.matrix}{the p x p matrix of correlation coefficients.}
+ \item{d}{values of d from the OU process for each trait.}
+ \item{B}{estimates of the regression coefficients, including intercepts. Coefficients are named according to the list U. For example, B1.2 is the coefficient corresponding to U[[1]][, 2], and if column 2 in U[[1]] is named "colname2", then the coefficient will be B1.colname2. Intercepts have the form B1.0.}
+ \item{B.se}{standard errors of the regression coefficients.}
+ \item{B.cov}{covariance matrix for regression coefficients.}
+ \item{B.zscore}{Z scores for the regression coefficients.}
+ \item{B.pvalue}{tests for the regression coefficients being different from zero.}
+ \item{logLik}{he log likelihood for either the restricted likelihood (REML = TRUE) or the overall likelihood (REML = FALSE).}
+ \item{AIC}{AIC for either the restricted likelihood (REML = TRUE) or the overall likelihood (REML = FALSE).}
+ \item{BIC}{BIC for either the restricted likelihood (REML = TRUE) or the overall likelihood (REML = FALSE).}
+ \item{REML}{whether REML is used rather than ML (TRUE or FALSE).}
+ \item{constrain.d}{whether or not values of d were constrained to be between 0 and 1 (TRUE or FALSE).}
+ \item{XX}{values of X in vectorized form, with each trait X[, i] standardized to have mean zero and standard deviation one.}
+ \item{UU}{design matrix with values in UU corresponding to XX; each variable U[[i]][, j] is standardized to have mean zero and standard deviation one.}
+ \item{MM}{vector of measurement standard errors corresponding to XX, with the standard errors suitably standardized.}
+ \item{Vphy}{the phylogenetic covariance matrix computed from phy and standardized to have determinant equal to one.}
+ \item{R}{covariance matrix of trait values relative to the standardized values of XX.}
+ \item{V}{overall estimated covariance matrix of residuals for XX including trait correlations, phylogenetic signal, and measurement error variances. This matrix can be used to simulate data for parametric bootstrapping. See examples.}
+ \item{C}{matrix V excluding measurement error variances.}
+ \item{convcode}{he convergence code provided by optim().}
+ \item{niter}{number of iterations performed by optim().}
+}
+\author{Anthony R. Ives}
+\references{
+ Zheng, L., A. R. Ives, T. Garland, B. R. Larget, Y. Yu, and K. F. Cao. 2009. New multivariate tests for phylogenetic signal and trait correlations applied to ecophysiological phenotypes of nine \emph{Manglietia} species. \emph{Functional Ecology} \bold{23}:1059--1069.
+
+}
+\examples{
+## Simple example using data without correlations or phylogenetic
+## signal. This illustrates the structure of the input data.
+
+phy <- rcoal(10, tip.label = 1:10)
+X <- matrix(rnorm(20), nrow = 10, ncol = 2)
+rownames(X) <- phy$tip.label
+U <- list(NULL, matrix(rnorm(10, mean = 10, sd = 4), nrow = 10, ncol = 1))
+rownames(U[[2]]) <- phy$tip.label
+SeM <- matrix(c(0.2, 0.4), nrow = 10, ncol = 2)
+rownames(SeM) <- phy$tip.label
+
+corphylo(X = X, SeM = SeM, U = U, phy = phy, method = "Nelder-Mead")
+
+\dontrun{
+## Simulation example for the correlation between two variables. The
+## example compares the estimates of the correlation coefficients from
+## corphylo when measurement error is incorporated into the analyses with
+## three other cases: (i) when measurement error is excluded, (ii) when
+## phylogenetic signal is ignored (assuming a "star" phylogeny), and (iii)
+## neither measurement error nor phylogenetic signal are included.
+
+## In the simulations, variable 2 is associated with a single
+## independent variable. This requires setting up a list U that has 2
+## elements: element U[[1]] is NULL and element U[[2]] is a n x 1 vector
+## containing simulated values of the independent variable.
+
+# Set up parameter values for simulating data
+n <- 50
+phy <- rcoal(n, tip.label = 1:n)
+
+R <- matrix(c(1, 0.7, 0.7, 1), nrow = 2, ncol = 2)
+d <- c(0.3, .95)
+B2 <- 1
+
+Se <- c(0.2, 1)
+SeM <- matrix(Se, nrow = n, ncol = 2, byrow = T)
+rownames(SeM) <- phy$tip.label
+
+# Set up needed matrices for the simulations
+p <- length(d)
+
+star <- stree(n)
+star$edge.length <- array(1, dim = c(n, 1))
+star$tip.label <- phy$tip.label
+
+Vphy <- vcv(phy)
+Vphy <- Vphy/max(Vphy)
+Vphy <- Vphy/exp(determinant(Vphy)$modulus[1]/n)
+
+tau <- matrix(1, nrow = n, ncol = 1) %*% diag(Vphy) - Vphy
+C <- matrix(0, nrow = p * n, ncol = p * n)
+for (i in 1:p) for (j in 1:p) {
+ Cd <- (d[i]^tau * (d[j]^t(tau)) * (1 - (d[i] * d[j])^Vphy))/(1 - d[i] * d[j])
+ C[(n * (i - 1) + 1):(i * n), (n * (j - 1) + 1):(j * n)] <- R[i, j] * Cd
+}
+MM <- matrix(SeM^2, ncol = 1)
+V <- C + diag(as.numeric(MM))
+
+## Perform a Cholesky decomposition of Vphy. This is used to generate
+## phylogenetic signal: a vector of independent normal random variables,
+## when multiplied by the transpose of the Cholesky deposition of Vphy will
+## have covariance matrix equal to Vphy.
+iD <- t(chol(V))
+
+# Perform Nrep simulations and collect the results
+Nrep <- 100
+cor.list <- matrix(0, nrow = Nrep, ncol = 1)
+cor.noM.list <- matrix(0, nrow = Nrep, ncol = 1)
+cor.noP.list <- matrix(0, nrow = Nrep, ncol = 1)
+cor.noMP.list <- matrix(0, nrow = Nrep, ncol = 1)
+d.list <- matrix(0, nrow = Nrep, ncol = 2)
+d.noM.list <- matrix(0, nrow = Nrep, ncol = 2)
+B.list <- matrix(0, nrow = Nrep, ncol = 3)
+B.noM.list <- matrix(0, nrow = Nrep, ncol = 3)
+B.noP.list <- matrix(0, nrow = Nrep, ncol = 3)
+for (rep in 1:Nrep) {
+ XX <- iD %*% rnorm(2 * n)
+ X <- matrix(XX, nrow = n, ncol = 2)
+ rownames(X) <- phy$tip.label
+
+ U <- list(NULL, matrix(rnorm(n, mean = 2, sd = 10), nrow = n, ncol = 1))
+ rownames(U[[2]]) <- phy$tip.label
+ colnames(U[[2]]) <- "V1"
+ X[,2] <- X[,2] + B2[1] * U[[2]][,1] - B2[1] * mean(U[[2]][,1])
+
+ z <- corphylo(X = X, SeM = SeM, U = U, phy = phy, method = "Nelder-Mead")
+ z.noM <- corphylo(X = X, U = U, phy = phy, method = "Nelder-Mead")
+ z.noP <- corphylo(X = X, SeM = SeM, U = U, phy = star, method = "Nelder-Mead")
+
+ cor.list[rep] <- z$cor.matrix[1, 2]
+ cor.noM.list[rep] <- z.noM$cor.matrix[1, 2]
+ cor.noP.list[rep] <- z.noP$cor.matrix[1, 2]
+ cor.noMP.list[rep] <- cor(cbind(lm(X[,1] ~ 1)$residuals, lm(X[,2] ~ U[[2]])$residuals))[1,2]
+
+ d.list[rep, ] <- z$d
+ d.noM.list[rep, ] <- z.noM$d
+
+ B.list[rep, ] <- z$B
+ B.noM.list[rep, ] <- z.noM$B
+ B.noP.list[rep, ] <- z.noP$B
+
+ show(c(rep, z$convcode, z$cor.matrix[1, 2], z$d))
+}
+correlation <- rbind(R[1, 2], mean(cor.list), mean(cor.noM.list),
+ mean(cor.noP.list), mean(cor.noMP.list))
+rownames(correlation) <- c("True", "With SeM and Phy", "Without SeM",
+ "Without Phy", "Without Phy or SeM")
+correlation
+
+signal.d <- rbind(d, colMeans(d.list), colMeans(d.noM.list))
+rownames(signal.d) <- c("True", "With SeM and Phy", "Without SeM")
+signal.d
+
+est.B <- rbind(c(0, 0, B2), colMeans(B.list), colMeans(B.noM.list),
+ colMeans(B.noP.list))
+rownames(est.B) <- c("True", "With SeM and Phy", "Without SeM", "Without Phy")
+colnames(est.B) <- rownames(z$B)
+est.B
+
+# Example simulation output
+# correlation
+ # [,1]
+# True 0.7000000
+# With SeM and Phy 0.7055958
+# Without SeM 0.3125253
+# Without Phy 0.4054043
+# Without Phy or SeM 0.3476589
+
+# signal.d
+ # [,1] [,2]
+# True 0.300000 0.9500000
+# With SeM and Phy 0.301513 0.9276663
+# Without SeM 0.241319 0.4872675
+
+# est.B
+ # B1.0 B2.0 B2.V1
+# True 0.00000000 0.0000000 1.0000000
+# With SeM and Phy -0.01285834 0.2807215 0.9963163
+# Without SeM 0.01406953 0.3059110 0.9977796
+# Without Phy 0.02139281 0.3165731 0.9942140
+}}
+\keyword{regression}
diff --git a/man/def.Rd b/man/def.Rd
index e47cd76..3b51382 100644
--- a/man/def.Rd
+++ b/man/def.Rd
@@ -39,7 +39,7 @@ given with the \code{\dots} (either "black" or 1).
If \code{regexp = TRUE} is used, then the names of the statements must be
quoted, e.g.:
-\code{def(tr$tip.label, "^Pan_" = "red")}
+\code{def(tr$tip.label, "^Pan_" = "red", regexp = TRUE)}
will return "red" for all labels starting with "Pan_".
}
diff --git a/man/del.gaps.Rd b/man/del.gaps.Rd
index 7a9ff18..59039b0 100644
--- a/man/del.gaps.Rd
+++ b/man/del.gaps.Rd
@@ -1,27 +1,33 @@
\name{del.gaps}
\alias{del.gaps}
-\title{
- Delete Alignment Gaps in DNA Sequences
+\alias{del.colgapsonly}
+\title{Delete Alignment Gaps in DNA Sequences}
+\description{
+ These functions remove indel gaps (\code{"-"}) in a sample of DNA
+ sequences.
}
\usage{
del.gaps(x)
+del.colgapsonly(x)
}
\arguments{
- \item{x}{a matrix, a list, or a vector containing the DNA sequences.}
-}
-\description{
- This function removes the insertion gaps (\code{"-"}) in a sample of
- DNA sequences.
+ \item{x}{a matrix, a list, or a vector containing the DNA
+ sequences; only matrices for \code{del.colgapsonly}.}
}
\details{
- The sequences can be either in \code{"DNAbin"} or in character format,
- but the returned object is always of class \code{"DNAbin"}. If
- \code{x} is a vector, then a vector is returned; if it is a list or a
- matrix, then a list is returned.
+ \code{del.gaps} remove all gaps, so the returned sequences may not
+ have all the same lengths and are therefore returned in a list.
+
+ \code{del.colgapsonly} removes the columns that contain only gaps
+ (useful when a small matrix is extracted from a large alignment).
+
+ The sequences can be either in \code{"DNAbin"} or in another format,
+ but the returned object is always of class \code{"DNAbin"}.
}
\value{
- A vector (if there is only one input sequence) or a list of class
- \code{"DNAbin"}.
+ \code{del.gaps} returns a vector (if there is only one input sequence)
+ or a list of class \code{"DNAbin"}; \code{del.colgapsonly} returns a
+ matrix.
}
\author{Emmanuel Paradis}
\seealso{
diff --git a/man/dist.topo.Rd b/man/dist.topo.Rd
index d1542ad..6c7814c 100644
--- a/man/dist.topo.Rd
+++ b/man/dist.topo.Rd
@@ -1,19 +1,21 @@
\name{dist.topo}
\alias{dist.topo}
\title{Topological Distances Between Two Trees}
+\description{
+ This function computes the topological distance between two
+ phylogenetic trees or among trees in a list (if \code{y = NULL} using
+ different methods.
+}
\usage{
-dist.topo(x, y, method = "PH85")
+dist.topo(x, y = NULL, method = "PH85")
}
\arguments{
- \item{x}{an object of class \code{"phylo"}.}
- \item{y}{an object of class \code{"phylo"}.}
+ \item{x}{an object of class \code{"phylo"} or of class
+ \code{"multiPhylo"}.}
+ \item{y}{an (optional) object of class \code{"phylo"}.}
\item{method}{a character string giving the method to be used: either
\code{"PH85"}, or \code{"score"}.}
}
-\description{
- This function computes the topological distance between two
- phylogenetic trees using different methods.
-}
\value{
a single numeric value.
}
@@ -65,6 +67,6 @@ dist.topo(x, y, method = "PH85")
ta <- rtree(30)
tb <- rtree(30)
dist.topo(ta, ta) # = 0
-dist.topo(ta, tb) # This is unlikely to be 0 !
+dist.topo(ta, tb) # This is unlikely to be 0!
}
\keyword{manip}
diff --git a/man/is.binary.tree.Rd b/man/is.binary.tree.Rd
index 992091b..1463da4 100644
--- a/man/is.binary.tree.Rd
+++ b/man/is.binary.tree.Rd
@@ -25,7 +25,7 @@ is.binary.tree(phy)
\seealso{
\code{\link{is.rooted}}, \code{\link{is.ultrametric}}, \code{\link{multi2di}}
}
-\author{Korbinian Strimmer (\url{http://www.stat.uni-muenchen.de/~strimmer/})}
+\author{Korbinian Strimmer}
\examples{
data("hivtree.newick") # example tree in NH format
tree.hiv <- read.tree(text = hivtree.newick)
diff --git a/man/landplants.Rd b/man/landplants.Rd
index 837a274..2b870e3 100644
--- a/man/landplants.Rd
+++ b/man/landplants.Rd
@@ -14,8 +14,7 @@ data(landplants.newick)
}
\source{
This tree is described in Sanderson (1997) and is also a
- data example in the software package r8s
- (\url{http://ginger.ucdavis.edu/r8s/}).
+ data example in the software package r8s.
}
\references{
Sanderson, M. J. (1997) A nonparametric approach to estimating
diff --git a/man/node.depth.Rd b/man/node.depth.Rd
index 925f9e1..73d6b97 100644
--- a/man/node.depth.Rd
+++ b/man/node.depth.Rd
@@ -2,7 +2,6 @@
\alias{node.depth}
\alias{node.depth.edgelength}
\alias{node.height}
-\alias{node.height.clado}
\title{Depth and Heights of Nodes and Tips}
\description{
These functions return the depths or heights of nodes and tips.
@@ -11,7 +10,6 @@
node.depth(phy, method = 1)
node.depth.edgelength(phy)
node.height(phy, clado.style = FALSE)
-node.height.clado(phy)
}
\arguments{
\item{phy}{an object of class "phylo".}
@@ -30,17 +28,11 @@ node.height.clado(phy)
\code{node.height} computes the heights of nodes and tips as plotted
by a phylogram or a cladogram.
-
- \code{node.height.clado} does the same but for a cladogram (this
- function will be removed soon as \code{node.height(phy, clado.style =
- TRUE)} does the same thing).
}
\value{
A numeric vector indexed with the node numbers of the matrix `edge' of
\code{phy}.
}
\author{Emmanuel Paradis}
-\seealso{
- \code{\link{plot.phylo}}
-}
+\seealso{\code{\link{plot.phylo}}}
\keyword{manip}
diff --git a/man/read.dna.Rd b/man/read.dna.Rd
index ca0a5c0..caba0d2 100644
--- a/man/read.dna.Rd
+++ b/man/read.dna.Rd
@@ -85,11 +85,7 @@ read.FASTA(file)
\code{read.FASTA} always returns a list of class \code{"DNAbin"}.
}
\references{
- Anonymous. FASTA format description.
- \url{http://www.ncbi.nlm.nih.gov/BLAST/fasta.html}
-
- Anonymous. IUPAC ambiguity codes.
- \url{http://www.ncbi.nlm.nih.gov/SNP/iupac.html}
+ Anonymous. FASTA format. \url{http://en.wikipedia.org/wiki/FASTA_format}
Felsenstein, J. (1993) Phylip (Phylogeny Inference Package) version
3.5c. Department of Genetics, University of Washington.
diff --git a/man/read.tree.Rd b/man/read.tree.Rd
index 4ad7e91..0ce1646 100644
--- a/man/read.tree.Rd
+++ b/man/read.tree.Rd
@@ -82,7 +82,7 @@ read.tree(file = "", text = NULL, tree.names = NULL, skip = 0,
\url{http://evolution.genetics.washington.edu/phylip/newick_doc.html}
Paradis, E. (2008) Definition of Formats for Coding Phylogenetic Trees
- in R. \url{http://ape.mpl.ird.fr/misc/FormatTreeR_28July2008.pdf}
+ in R. \url{http://ape-package.ird.fr/misc/FormatTreeR_24Oct2012.pdf}
Paradis, E. (2012) \emph{Analysis of Phylogenetics and Evolution with
R (Second Edition).} New York: Springer.
diff --git a/man/rlineage.Rd b/man/rlineage.Rd
index e01a810..1f7299f 100644
--- a/man/rlineage.Rd
+++ b/man/rlineage.Rd
@@ -1,13 +1,16 @@
\name{rlineage}
\alias{rlineage}
\alias{rbdtree}
+\alias{rphylo}
\alias{drop.fossil}
\title{Tree Simulation Under the Time-Dependent Birth--Death Models}
\description{
- These two functions simulate phylogenies under any time-dependent
- birth--death model. \code{lineage} generates a complete tree including
- the species that go extinct; \code{rbdtree} generates a tree with only
- the species until present; \code{drop.fossil} is a utility function to
+ These three functions simulate phylogenies under any time-dependent
+ birth--death model: \code{rlineage} generates a complete tree including
+ the species going extinct before present; \code{rbdtree} generates a
+ tree with only the species living at present (thus the tree is
+ ultrametric); \code{rphylo} generates a tree with a fixed number of
+ species at present time. \code{drop.fossil} is a utility function to
remove the extinct species.
}
\usage{
@@ -15,6 +18,8 @@ rlineage(birth, death, Tmax = 50, BIRTH = NULL,
DEATH = NULL, eps = 1e-6)
rbdtree(birth, death, Tmax = 50, BIRTH = NULL,
DEATH = NULL, eps = 1e-6)
+rphylo(n, birth, death, BIRTH = NULL, DEATH = NULL,
+ T0 = 50, fossils = FALSE, eps = 1e-06)
drop.fossil(phy, tol = 1e-8)
}
\arguments{
@@ -23,23 +28,39 @@ drop.fossil(phy, tol = 1e-8)
\item{Tmax}{a numeric value giving the length of the simulation.}
\item{BIRTH, DEATH}{a (vectorized) function which is the primitive
of \code{birth} or \code{death}. This can be used to speed-up the
- computation. By default, numerical integration is done.}
+ computation. By default, a numerical integration is done.}
\item{eps}{a numeric value giving the time resolution of the
simulation; this may be increased (e.g., 0.001) to shorten
computation times.}
+ \item{n}{the number of species living at present time.}
+ \item{T0}{the time at present (for the backward-in-time algorithm).}
+ \item{fossils}{a logical value specifying whether to output the
+ lineages going extinct.}
\item{phy}{an object of class \code{"phylo"}.}
\item{tol}{a numeric value giving the tolerance to consider a species
as extinct.}
}
\details{
- Both functions use continuous-time algorithms described in the
- references. The models are time-dependent birth--death models as
- described in Kendall (1948). Speciation (birth) and extinction (death)
- rates may be constant or vary through time according to an \R function
- specified by the user. In the latter case, \code{BIRTH} and/or
- \code{DEATH} may be used if the primitives of \code{birth} and
- \code{death} are known. In these functions time is the formal argument
- and must be named \code{t}.
+ These three functions use continuous-time algorithms: \code{rlineage}
+ and \code{rbdtree} use the forward-in-time algorithms described in
+ Paradis (2011), whereas \code{rphylo} uses a backward-in-time
+ algorithm from Stadler (2011). The models are time-dependent
+ birth--death models as described in Kendall (1948). Speciation
+ (birth) and extinction (death) rates may be constant or vary through
+ time according to an \R function specified by the user. In the latter
+ case, \code{BIRTH} and/or \code{DEATH} may be used if the primitives
+ of \code{birth} and \code{death} are known. In these functions time is
+ the formal argument and must be named \code{t}.
+
+ Note that \code{rphylo} simulates trees in a way similar to what
+ the package \pkg{TreeSim} does, the difference is in the
+ parameterization of the time-dependent models which is here the same
+ than used in the two other functions. In this parameterization scheme,
+ time is measured from past to present (see details in Paradis 2015
+ which includes a comparison of these algorithms).
+
+ The difference between \code{rphylo} and \code{rphylo(... fossils
+ = TRUE)} is the same than between \code{rbdtree} and \code{rlineage}.
}
\value{
An object of class \code{"phylo"}.
@@ -51,6 +72,12 @@ drop.fossil(phy, tol = 1e-8)
Paradis, E. (2011) Time-dependent speciation and extinction from
phylogenies: a least squares approach. \emph{Evolution}, \bold{65},
661--672.
+
+ Paradis, E. (2015) Random phylogenies and the distribution of
+ branching times. Manuscript.
+
+ Stadler, T. (2011) Simulating trees with a fixed number of extant
+ species. \emph{Systematic Biology}, \bold{60}, 676--684.
}
\author{Emmanuel Paradis}
\seealso{
diff --git a/man/root.Rd b/man/root.Rd
index 05ef593..499605c 100644
--- a/man/root.Rd
+++ b/man/root.Rd
@@ -51,6 +51,17 @@ is.rooted(phy)
root, or if there is a \code{root.edge} element. In all other cases,
\code{is.rooted} returns \code{FALSE}.
}
+\note{
+ The use of \code{resolve.root = TRUE} together with \code{node = }
+ gives an error if the specified node is the current root of the
+ tree. This is because there is an ambiguity when resolving a node in
+ an unrooted tree if no explicit outgroup. If the node is not the
+ current root, the ambiguity is solved arbitrarily by considering the
+ clade on the right of \code{node} (when the tree is plotted by
+ default) as the ingroup. See a detailed explanation there:
+
+ \url{https://stat.ethz.ch/pipermail/r-sig-phylo/2015-April/003979.html}.
+}
\value{
an object of class \code{"phylo"} for \code{root} and \code{unroot}; a
single logical value for \code{is.rooted}.
diff --git a/man/seg.sites.Rd b/man/seg.sites.Rd
index 81f046a..5ef6c43 100644
--- a/man/seg.sites.Rd
+++ b/man/seg.sites.Rd
@@ -22,11 +22,6 @@ seg.sites(x)
sites.
}
\author{Emmanuel Paradis}
-\note{
- The present version looks for the sites which are ``variable'' in the
- data in terms of different \emph{letters}. This may give unexpected
- results if there are ambiguous bases in the data.
-}
\seealso{
\code{\link{base.freq}}, \code{\link[pegas]{theta.s}},
\code{\link[pegas]{nuc.div}}
diff --git a/man/skyline.Rd b/man/skyline.Rd
index 718ff22..50d373e 100644
--- a/man/skyline.Rd
+++ b/man/skyline.Rd
@@ -4,7 +4,6 @@
\alias{skyline.coalescentIntervals}
\alias{skyline.collapsedIntervals}
\alias{find.skyline.epsilon}
-
\title{Skyline Plot Estimate of Effective Population Size}
\usage{
skyline(x, \dots)
@@ -21,35 +20,27 @@ find.skyline.epsilon(ci, GRID=1000, MINEPS=1e-6, \dots)
\item{epsilon}{collapsing parameter that controls the amount of smoothing
(allowed range: from \code{0} to \code{ci$total.depth}, default value: 0). This is the same parameter as in
\link{collapsed.intervals}.}
-
\item{old.style}{Parameter to choose between two slightly different variants of the
generalized skyline plot (Strimmer and Pybus, pers. comm.). The default value \code{FALSE} is
recommended.}
-
\item{ci}{coalescent intervals (i.e. an object of class \code{"coalescentIntervals"})}
-
\item{GRID}{Parameter for the grid search for \code{epsilon} in \code{find.skyline.epsilon}.}
-
\item{MINEPS}{Parameter for the grid search for \code{epsilon} in \code{find.skyline.epsilon}.}
-
\item{\dots}{Any of the above parameters.}
-
}
\description{
-
\code{skyline} computes the \emph{generalized skyline plot} estimate of effective population size
- from an estimated phylogeny. The demographic history is approximated by
+ from an estimated phylogeny. The demographic history is approximated by
a step-function. The number of parameters of the skyline plot (i.e. its smoothness)
- is controlled by a parameter \code{epsilon}.
-
+ is controlled by a parameter \code{epsilon}.
+
\code{find.skyline.epsilon} searches for an optimal value of the \code{epsilon} parameter,
i.e. the value that maximizes the AICc-corrected log-likelihood (\code{logL.AICc}).
}
-
\details{
-\code{skyline} implements the \emph{generalized skyline plot} introduced in
+\code{skyline} implements the \emph{generalized skyline plot} introduced in
Strimmer and Pybus (2001). For \code{epsilon = 0} the
-generalized skyline plot degenerates to the
+generalized skyline plot degenerates to the
\emph{classic skyline plot} described in
Pybus et al. (2000). The latter is in turn directly related to lineage-through-time plots
(Nee et al., 1995).
@@ -61,31 +52,21 @@ Pybus et al. (2000). The latter is in turn directly related to lineage-through-
\item{time}{ A vector with the time at the end of each coalescent
interval (i.e. the accumulated interval lengths from the beginning of the first interval
to the end of an interval)}
-
- \item{interval.length}{ A vector with the length of each
- interval.}
-
+ \item{interval.length}{ A vector with the length of each interval.}
\item{population.size}{A vector with the effective population size of each interval.}
-
- \item{parameter.count}{ Number of free parameters in the skyline plot.}
+ \item{parameter.count}{ Number of free parameters in the skyline plot.}
\item{epsilon}{The value of the underlying smoothing parameter.}
-
\item{logL}{Log-likelihood of skyline plot (see Strimmer and Pybus, 2001).}
-
\item{logL.AICc}{AICc corrected log-likelihood (see Strimmer and Pybus, 2001).}
-\code{find.skyline.epsilon} returns the value of the \code{epsilon} parameter
- that maximizes \code{logL.AICc}.
+ \code{find.skyline.epsilon} returns the value of the \code{epsilon} parameter
+ that maximizes \code{logL.AICc}.
}
-
-\author{Korbinian Strimmer (\url{http://www.stat.uni-muenchen.de/~strimmer/})}
-
+\author{Korbinian Strimmer}
\seealso{
\code{\link{coalescent.intervals}}, \code{\link{collapsed.intervals}},
\code{\link{skylineplot}}, \code{\link{ltt.plot}}.
}
-
-
\references{
Strimmer, K. and Pybus, O. G. (2001) Exploring the demographic history
of DNA sequences using the generalized skyline plot. \emph{Molecular
@@ -100,7 +81,6 @@ Pybus et al. (2000). The latter is in turn directly related to lineage-through-
Transactions of the Royal Society of London. Series B. Biological
Sciences}, \bold{349}, 25--31.
}
-
\examples{
# get tree
data("hivtree.newick") # example tree in NH format
@@ -113,7 +93,6 @@ ci <- coalescent.intervals(tree.hiv) # from tree
cl1 <- collapsed.intervals(ci,0)
cl2 <- collapsed.intervals(ci,0.0119)
-
#### classic skyline plot ####
sk1 <- skyline(cl1) # from collapsed intervals
sk1 <- skyline(ci) # from coalescent intervals
@@ -134,7 +113,6 @@ sk2
plot(sk2)
-
# classic and generalized skyline plot together in one plot
plot(sk1, show.years=TRUE, subst.rate=0.0023, present.year = 1997, col=c(grey(.8),1))
lines(sk2, show.years=TRUE, subst.rate=0.0023, present.year = 1997)
diff --git a/man/skylineplot.Rd b/man/skylineplot.Rd
index f9e7166..1f9db7e 100644
--- a/man/skylineplot.Rd
+++ b/man/skylineplot.Rd
@@ -3,7 +3,6 @@
\alias{plot.skyline}
\alias{lines.skyline}
\alias{skylineplot.deluxe}
-
\title{Drawing Skyline Plot Graphs}
\usage{
\method{plot}{skyline}(x, show.years=FALSE, subst.rate, present.year, \dots)
@@ -13,28 +12,19 @@ skylineplot.deluxe(tree, \dots)
}
\arguments{
\item{x}{skyline plot data (i.e. an object of class \code{"skyline"}).}
-
- \item{z}{Either an ultrametric tree (i.e. an object of class \code{"phylo"}),
+ \item{z}{Either an ultrametric tree (i.e. an object of class \code{"phylo"}),
or coalescent intervals (i.e. an object of class \code{"coalescentIntervals"}), or
collapsed coalescent intervals (i.e. an object of class \code{"collapsedIntervals"}).}
-
-
\item{tree}{ultrametric tree (i.e. an object of class \code{"phylo"}).}
-
-
\item{show.years}{option that determines whether the time is plotted in units of
of substitutions (default) or in years (requires specification of substution rate
and year of present).}
+ \item{subst.rate}{substitution rate (see option show.years).}
+ \item{present.year}{present year (see option show.years).}
-
- \item{subst.rate}{substitution rate (see option show.years).}
- \item{present.year}{present year (see option show.years).}
-
- \item{\dots}{further arguments to be passed on to \code{skyline()} and \code{plot()}.}
-
+ \item{\dots}{further arguments to be passed on to \code{skyline()} and \code{plot()}.}
}
\description{
-
These functions provide various ways to draw \emph{skyline plot} graphs
on the current graphical device. Note that \code{skylineplot(z, \dots)} is simply
a shortcut for \code{plot(skyline(z, \dots))}.
@@ -45,40 +35,32 @@ skylineplot.deluxe(tree, \dots)
\details{
See \code{\link{skyline}} for more details (incl. references) about the skyline plot method.
}
-
-
-\author{Korbinian Strimmer (\url{http://www.stat.uni-muenchen.de/~strimmer/})}
-
+\author{Korbinian Strimmer}
\seealso{
\code{\link[graphics]{plot}} and \code{\link[graphics]{lines}} for the basic plotting
function in R, \code{\link{coalescent.intervals}}, \code{\link{skyline}}
}
-
\examples{
# get tree
data("hivtree.newick") # example tree in NH format
tree.hiv <- read.tree(text = hivtree.newick) # load tree
-
#### classic skyline plot
skylineplot(tree.hiv) # shortcut
-
#### plot classic and generalized skyline plots and estimate epsilon
-sk.opt <- skylineplot.deluxe(tree.hiv)
+sk.opt <- skylineplot.deluxe(tree.hiv)
sk.opt$epsilon
-
#### classic and generalized skyline plot ####
-sk1 <- skyline(tree.hiv)
-sk2 <- skyline(tree.hiv, 0.0119)
+sk1 <- skyline(tree.hiv)
+sk2 <- skyline(tree.hiv, 0.0119)
# use years rather than substitutions as unit for the time axis
plot(sk1, show.years=TRUE, subst.rate=0.0023, present.year = 1997, col=c(grey(.8),1))
lines(sk2, show.years=TRUE, subst.rate=0.0023, present.year = 1997)
legend(.15,500, c("classic", "generalized"), col=c(grey(.8),1),lty=1)
-
#### various skyline plots for different epsilons
layout(mat= matrix(1:6,2,3,byrow=TRUE))
ci <- coalescent.intervals(tree.hiv)
diff --git a/man/where.Rd b/man/where.Rd
index 09e0ab8..1fd9402 100644
--- a/man/where.Rd
+++ b/man/where.Rd
@@ -18,8 +18,7 @@ where(x, pattern)
each sequence.
Patterns may be overlapping. For instance, if \code{pattern = "tata"}
- and the sequence starts with `tatata', then the vector returned will
- be c(1, 3).
+ and the sequence starts with `tatata', then the output will be c(1, 3).
}
\value{
a vector of integers or a list of such vectors.
diff --git a/man/woodmouse.Rd b/man/woodmouse.Rd
index d90ce66..dc92570 100644
--- a/man/woodmouse.Rd
+++ b/man/woodmouse.Rd
@@ -12,8 +12,7 @@
data(woodmouse)
}
\format{
- The data are stored in an ASCII file using the sequential format for
- DNA sequences which is read with `read.dna()'.
+ An object of class \code{"DNAbin"}.
}
\source{
Michaux, J. R., Magnanou, E., Paradis, E., Nieberding, C. and Libois,
diff --git a/man/write.dna.Rd b/man/write.dna.Rd
index 2ff245d..3d00a21 100644
--- a/man/write.dna.Rd
+++ b/man/write.dna.Rd
@@ -76,11 +76,7 @@ write.dna(x, file, format = "interleaved", append = FALSE,
}
\author{Emmanuel Paradis}
\references{
- Anonymous. FASTA format description.
- \url{http://www.ncbi.nlm.nih.gov/BLAST/fasta.html}
-
- Anonymous. IUPAC ambiguity codes.
- \url{http://www.ncbi.nlm.nih.gov/SNP/iupac.html}
+ Anonymous. FASTA format. \url{http://en.wikipedia.org/wiki/FASTA_format}
Felsenstein, J. (1993) Phylip (Phylogeny Inference Package) version
3.5c. Department of Genetics, University of Washington.
diff --git a/src/dist_dna.c b/src/dist_dna.c
index 8adcb68..650e02f 100644
--- a/src/dist_dna.c
+++ b/src/dist_dna.c
@@ -1,6 +1,6 @@
-/* dist_dna.c 2013-08-20 */
+/* dist_dna.c 2015-02-25 */
-/* Copyright 2005-2013 Emmanuel Paradis */
+/* Copyright 2005-2015 Emmanuel Paradis */
/* This file is part of the R-package `ape'. */
/* See the file ../COPYING for licensing issues. */
@@ -1071,20 +1071,18 @@ void BaseProportion(unsigned char *x, int *n, double *BF)
void SegSites(unsigned char *x, int *n, int *s, int *seg)
{
- int i, j;
- unsigned char basis;
+ int i, ib, j;
for (j = 0; j < *s; j++) {
- i = *n * j;
- while (!KnownBase(x[i])) i++;
- basis = x[i];
- i++;
- while (i < *n * (j + 1)) {
- if (!KnownBase(x[i]) || x[i] == basis) i++;
- else {
- seg[j] = 1;
- break;
+ for (i = *n * j; i < *n * (j + 1) - 1; i++) {
+ if (!KnownBase(x[i])) continue;
+ for (ib = i + 1; ib < *n * (j + 1); ib++) {
+ if (DifferentBase(x[i], x[ib])) {
+ seg[j] = 1;
+ break;
+ }
}
+ if (seg[j]) break;
}
}
}
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/r-cran-ape.git
More information about the debian-med-commit
mailing list