[med-svn] [SCM] vegan branch, master, updated. debian/2.0-3-1-10-g049dcc3

Charles Plessy plessy at debian.org
Thu Apr 4 05:44:28 UTC 2013


The following commit has been merged in the master branch:
commit ab1ff61d10270c13654a42bdd1d4fbf8df4e6452
Author: Charles Plessy <plessy at debian.org>
Date:   Thu Apr 4 14:41:52 2013 +0900

    Imported Upstream version 2.0-5

diff --git a/DESCRIPTION b/DESCRIPTION
index 05aa2dd..f16939f 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: vegan
 Title: Community Ecology Package
-Version: 2.0-4
-Date: June 18, 2012
+Version: 2.0-5
+Date: October 8, 2012
 Author: Jari Oksanen, F. Guillaume Blanchet, Roeland Kindt, Pierre
         Legendre, Peter R. Minchin, R. B. O'Hara, Gavin L. Simpson,
         Peter Solymos, M. Henry H. Stevens, Helene Wagner
@@ -13,6 +13,6 @@ Description: Ordination methods, diversity analysis and other functions
         for community and vegetation ecologists.
 License: GPL-2
 URL: http://cran.r-project.org, http://vegan.r-forge.r-project.org/
-Packaged: 2012-06-18 10:55:31 UTC; jarioksa
+Packaged: 2012-10-08 13:52:59 UTC; jarioksa
 Repository: CRAN
-Date/Publication: 2012-06-19 09:32:41
+Date/Publication: 2012-10-08 14:38:53
diff --git a/MD5 b/MD5
index 73a69e4..5439de7 100644
--- a/MD5
+++ b/MD5
@@ -1,6 +1,6 @@
-12748bc3d050f5cc3e9705f885fec6dd *DESCRIPTION
-b06b62a7c78c5d4ec84b2f567fd8e465 *NAMESPACE
-0f88ae5115f589efd416613a823cbb3f *R/AIC.radfit.R
+841d99790c24b6b46b9be2203bdc45fc *DESCRIPTION
+439a7709754560d1bdd238e351a324c0 *NAMESPACE
+4b8531b446af54510e5fb31f841aed2f *R/AIC.radfit.R
 5c5fdbcdc2a38e2cbafdb8f2c5eb2e08 *R/CCorA.R
 6592cf7dc692f87b4a147eb625e18624 *R/MDSrotate.R
 1a95c5873b1546683487e17aae4fe511 *R/MOStest.R
@@ -13,16 +13,16 @@ d80688d78aba3cd9367ffaaaec6ec252 *R/TukeyHSD.betadisper.R
 a7e4a96ac08502dc4cd7023d0510e94f *R/add1.cca.R
 bc63dc46c564df2dea69d90f25e52ba3 *R/ade2vegancca.R
 3fea698281bc0b4c3a5ad26f4d44d0e2 *R/adipart.R
-3018034f2134aef536ee9a5ef6b50a06 *R/adipart.default.R
-a71174cb8825cf6ed74612bc7cd41fcc *R/adipart.formula.R
-e71c569aecb494fdf872bbebe65276f9 *R/adonis.R
+ca2eae68396cb5a16355751efb906557 *R/adipart.default.R
+05387ee9e552fcec123b4b922e837eaa *R/adipart.formula.R
+39805e701d7a10abe111056a5b8afcfe *R/adonis.R
 38726b7dad817378b3fc41094c11595a *R/alias.cca.R
 fa03b3668c0435eebc778b9def908be9 *R/anosim.R
 a4f23289c4a5eab2a3587292b306d497 *R/anova.betadisper.R
 a3767eb78bd5d62341ff198fbf65d1b1 *R/anova.cca.R
-70ca1744d8ba7f7601f8cccd826838d0 *R/anova.ccabyaxis.R
+363461c0bab2687ece5bf521398db928 *R/anova.ccabyaxis.R
 350cc1a6e930a195c80b5c7aef0d9b7e *R/anova.ccabymargin.R
-42491e7011d7ae79ea2a8db2086768f9 *R/anova.ccabyterm.R
+145b004abb6a5dca32fc8f60338af590 *R/anova.ccabyterm.R
 97cbe54d0f4f7adee7a20b6c982d1ecf *R/anova.ccanull.R
 7fab08bcc596df60a22c4b04c8507121 *R/anova.prc.R
 6fb2bf929aed44ef41bfd4dfc6e010cc *R/as.fisher.R
@@ -36,14 +36,14 @@ a7f01bd69394d5554cf10279a2690080 *R/as.preston.R
 124345e75a121a5d750d1b306156614a *R/as.ts.oecosimu.R
 251aa91b87e2f964f24fe61791eedac9 *R/as.ts.permat.R
 fbec6d133dea10372ce082c7035a8ab2 *R/beals.R
-36d76d3619d9effc5819e907a45d7476 *R/betadisper.R
+d25853e603f6af59234db7714c55f71f *R/betadisper.R
 92d22ad6b8870005ed86940cd90f2077 *R/betadiver.R
 46ae3f75a0b483fecab589637d72a307 *R/bgdispersal.R
 4603ea944d470a9e284cb6cab6d75529 *R/bioenv.R
 4018f8ddbc314125dbd17d44f710147e *R/bioenv.default.R
 82a17642aae8182c9345a03a846a853c *R/bioenv.formula.R
 875d40515bf55ee65dc7fcdefb9f52d1 *R/biplot.CCorA.R
-329ff7c5c96cbfe243b933a748af979c *R/biplot.rda.R
+e83522ded9481ebde69e61419d0033b7 *R/biplot.rda.R
 0999bb90f22b72fade2ca6adbd01758f *R/boxplot.betadisper.R
 f5abc1a3e5417e53a78cf2054f46d0a6 *R/boxplot.specaccum.R
 cbf54233db3c2839101f98e02eb538dd *R/bstick.R
@@ -52,7 +52,7 @@ cbf54233db3c2839101f98e02eb538dd *R/bstick.R
 b5a4874f7763f1793b3068eec4e859d5 *R/bstick.default.R
 4bf61eddbdd1cec0bd03150cd478f3d9 *R/bstick.prcomp.R
 2ad6f2071ad822c631e04705df2b245c *R/bstick.princomp.R
-2da3da8e0214c1a9a3843c47782695a4 *R/cIndexKM.R
+b3ec7384552e94c46102d26e38c40538 *R/cIndexKM.R
 a6df607186ceb18d204494b6a33816d4 *R/calibrate.R
 f093f401495379d30354e6424606190a *R/calibrate.cca.R
 f56b52d53b17c7dc8d8c9accd5a3401a *R/calibrate.ordisurf.R
@@ -65,7 +65,7 @@ e01e3acecdb9ac8d9195937e9879d126 *R/cca.formula.R
 c66d8fbe69ccca94f2ee8f777ff16ae2 *R/checkSelect.R
 6c516cbb81d76ad6c22175809211eec4 *R/clamtest.R
 365fa22c822e835218740e858d256526 *R/coef.cca.R
-3a7297a8e2b2858f3b776a35689e5eca *R/coef.radfit.R
+8ed92f141523fe0cc4c3bb51807b8b07 *R/coef.radfit.R
 ea10763445cb76b219d18bb274109df5 *R/coef.rda.R
 ab87ce0f23c88b6b40207a7597fa9b64 *R/commsimulator.R
 722959743928d23213409c778c6acbc2 *R/confint.MOStest.R
@@ -75,8 +75,9 @@ d0f10f160ac99ba936824a49c608868a *R/cophenetic.spantree.R
 edee3aaced61290b219985d0ce69155c *R/coverscale.R
 1b1a6343072d69c5ccbf9a72ba068cbd *R/decorana.R
 1169ef2aa4cd76c33c6166b1b253a665 *R/decostand.R
-28b0ea3a91ffd243e8da8d9430978cc9 *R/density.oecosimu.R
-49baf1d517af98a4d996c6a26886de48 *R/densityplot.oecosimu.R
+72e1985e74ccaf8785e7d6b6f2005018 *R/density.anosim.R
+4a13947927b175862e2266ff9589f2a0 *R/density.oecosimu.R
+b18236e376b3e0fb38faa1cf10139f30 *R/densityplot.oecosimu.R
 e5b54fa580331ab24d28dc59110c45fe *R/designdist.R
 8fb0105cb21a5e425f72e0085fa669a2 *R/deviance.cca.R
 52f41185396ddd9acdcee9df7298d65a *R/deviance.rda.R
@@ -105,16 +106,17 @@ d9ca8944dc5fdfc396e50d118c2638c3 *R/fitspecaccum.R
 fdd1e11bb4f5ed4c31497dd3dfe18880 *R/fitted.capscale.R
 1c356de7c7721d2a6cb04abf54d91274 *R/fitted.cca.R
 0080b65cfd48bac5e53961b8e12682e5 *R/fitted.procrustes.R
-5f97ed9b7659f0348db8ca29a119c1bf *R/fitted.radfit.R
+bdb5429c7c23d1f730c5c7c938fb5e09 *R/fitted.radfit.R
 6b3b45debf635860ae9c1f3bc17686fc *R/fitted.rda.R
 57c9a7ccff6a9c066b2aba3475f2330b *R/goodness.R
 c3bc6f9771503a4b93f4dce3096fd874 *R/goodness.cca.R
 7af5f06020065621d8466decb16e0aa4 *R/goodness.metaMDS.R
 3676ae4dd5c789334fb822a8caa81979 *R/goodness.rda.R
 8a767726c40223a58d4055759bf41efe *R/head.summary.cca.R
+d17f4f6be45b52e01cd605b09b56a80a *R/hierParseFormula.R
 3d19236ee5dd2f1c678061773895e86f *R/hiersimu.R
-14cebf89aebdc5c9b4c910dee6417487 *R/hiersimu.default.R
-9cc65f22f3489d18c27fc07df9f32337 *R/hiersimu.formula.R
+5c12aa3da8e9ed06db19bad2f1eea41a *R/hiersimu.default.R
+edf53c3358944421756412b991432bd7 *R/hiersimu.formula.R
 d02fc9c672a9b2c4a31065702a3381be *R/humpfit.R
 1637bd10b39801c14b65656f29dafcf1 *R/identify.ordiplot.R
 ba2579e418e11343a471548b2f2bc44a *R/indpower.R
@@ -125,17 +127,18 @@ c63972a171f76f92652feeb2daf30e82 *R/isomap.R
 1e167e69edcee4aa651d97bef81b31e9 *R/isomapdist.R
 5abdcd58cf3811e482543d5207114331 *R/kendall.global.R
 814a92ded8985a249d4f5973f2c1fd6c *R/kendall.post.R
+14f4282484b5ff091f36b30e190693f8 *R/labels.envfit.R
 131b9e9e84f94acdf044430143e19ac8 *R/lines.humpfit.R
 f1d30acca998f0fe17e8363203f1b920 *R/lines.permat.R
 eb4e11e71eeefa6ec64e4a2580b8af75 *R/lines.prestonfit.R
 a5d9b7fa31477efc0a1ff76a62988a8e *R/lines.procrustes.R
-2f7580776c503a92034de7a62a271d0c *R/lines.radline.R
+39604c069428cda7c9d2ed199ac4e28a *R/lines.radline.R
 f58c42977c7b03c9d2252295068b8846 *R/lines.spantree.R
 e0f782c5a0cb558189b5b15a5bea072f *R/linestack.R
 1dcc7e0504b5468a3bb2253924901e50 *R/make.cepnames.R
-409328452f8b9d78bc5a2ed4499559b8 *R/mantel.R
+106633a42f577ab3e266b146ba4447e0 *R/mantel.R
 fdb2f4786b31866197c80d827584edaf *R/mantel.correlog.R
-3109021e5461bac44124343acfd3ec81 *R/mantel.partial.R
+6ba992d75c8ec43717adf957c8ff2427 *R/mantel.partial.R
 e054f13ad65a7f2616561c73557b412b *R/meandist.R
 57cb748570098b7e5a5aedbddb39fb84 *R/metaMDS.R
 26b26e400ead4cf3de31d7eab29c6984 *R/metaMDSdist.R
@@ -148,8 +151,8 @@ beae6832197823e4ade3569c021b1693 *R/mrpp.R
 0c8ef224eeeb965c8b60bb63d5adf10e *R/mso.R
 0595a503a1542927040f1b643e85d9a2 *R/msoplot.R
 7c219818ce5841e957db47f16986080b *R/multipart.R
-569dd4c9ebbe2ab60d33ff3e4e4bda38 *R/multipart.default.R
-dde94f190e182c04ffa012e3a2ba38cd *R/multipart.formula.R
+1757260e533d3b5521c86099be28dc5a *R/multipart.default.R
+4f3e2c82d5783c04f9a50761c82e2f02 *R/multipart.formula.R
 f5e79cb1c2dc1fcabb6e6b5cb4dc0828 *R/nestedbetasor.R
 85d4744650c1e2a0edf16809b77f7ab4 *R/nestedchecker.R
 c15884dd28790c7521ecb33175c86e5c *R/nesteddisc.R
@@ -158,7 +161,7 @@ e65023174f4ce8874a2f88f332af5a37 *R/nestedn0.R
 e27188404c8386162d6af0fb47110a84 *R/nestedtemp.R
 74b2723851155de631716fa479f8ea38 *R/no.shared.R
 47973ff187f68836a19d20ea37c60868 *R/nobs.R
-d5e823fadab8ffd703b7e8d0c7d3232b *R/oecosimu.R
+20078fe8d961f3da96a0a589fcd53ee8 *R/oecosimu.R
 7b3988a207ecfe1ea574c5857ffcd2a3 *R/orderingKM.R
 e3d108eed97633040fa22c2b384e19e4 *R/ordiArgAbsorber.R
 d2d5d94504676e6b1edf25b0c024edf3 *R/ordiArrowMul.R
@@ -200,7 +203,7 @@ c53e9402a842833d80a8df39c0adee6f *R/orglpoints.R
 e08110689dfeb1098cb4d9194f084c66 *R/permatfull.R
 26a9634c5ad6bc16e2e24c283e33b761 *R/permatswap.R
 909306255cee4f36d2ba7ba13d376e90 *R/permuted.index.R
-2800b02ade703d3fdcff1a86380c53e7 *R/permutest.betadisper.R
+7aedada06df1c5e0dff74f77e4479fbb *R/permutest.betadisper.R
 4764a3d49455270e5217b72aa4d68787 *R/permutest.cca.R
 b4e77b98f86c4b567d687b64e3aa8812 *R/persp.renyiaccum.R
 011a26868189ef4f3516a1b1931a2ea1 *R/persp.tsallisaccum.R
@@ -213,7 +216,7 @@ d6efc684c3506870243a337fadd85ed8 *R/plot.cca.R
 38ccde16c9eb9028219d27f14d343b3e *R/plot.clamtest.R
 8c043a9b7262c33ec2054045cdaa1811 *R/plot.contribdiv.R
 0ab3b62ac155ede193867f43640dbd34 *R/plot.decorana.R
-0b0701c8bdeb6695a75503b6c0e93466 *R/plot.envfit.R
+090911b19b4659e69ba80dab85bc7ff7 *R/plot.envfit.R
 10bf121300b684a8173f680de54f452a *R/plot.fisherfit.R
 9a4f1746e6e5b80b48994f404e72eb74 *R/plot.humpfit.R
 ed258eefbe3facce3533a16395217fab *R/plot.isomap.R
@@ -231,7 +234,7 @@ fdc1beae72f52a43883861a8b56bf289 *R/plot.prc.R
 11d54a7950f3e00548eb8c32f4f5be98 *R/plot.procrustes.R
 c7e5c6d58944d75ab6dac163e051769f *R/plot.profile.fisherfit.R
 02ff38f3fb337a63534356255b8641a9 *R/plot.rad.R
-9935e5d441eda4ec3445ae9f3eaa4812 *R/plot.radfit.R
+fc2dc1b63ae6f50067a7a376c736394b *R/plot.radfit.R
 97e04e44547a38580a988bbd684e4cb3 *R/plot.radfit.frame.R
 360dec911e8d4e772f888d89b8e0f6f7 *R/plot.radline.R
 8da525de886af7ced0065ecc8add6fa9 *R/plot.renyi.R
@@ -248,7 +251,7 @@ a0e1e2d579fa8c1992a26a2e8d435750 *R/points.metaMDS.R
 a54bcddf1b7a44ee1f86ae4eaccb7179 *R/points.ordiplot.R
 e352171f478eb27cf4a875cc3a1693fc *R/points.orditkplot.R
 7409704e2e94cd051524e8c5af3bdcb4 *R/points.procrustes.R
-d76776427a3f064c91b4c1c14da0d96b *R/points.radline.R
+80d9cee7ff1fa7ee8cb18850711a14b2 *R/points.radline.R
 b4fbbb0786258e1e83c4262e0db2aa43 *R/poolaccum.R
 91aa7fd2fbd99f8e325932d59886dac7 *R/postMDS.R
 f9dcd972e5c81ce936c9ec5b296d484c *R/prc.R
@@ -256,6 +259,7 @@ f9dcd972e5c81ce936c9ec5b296d484c *R/prc.R
 049f41cca1b39bf0a221723855cffcff *R/predict.decorana.R
 e9107f2a97c429ff8e45c62a3671c281 *R/predict.fitspecaccum.R
 06cca728e43d29da2528b01dccb26962 *R/predict.humpfit.R
+3eaaaf25580077e7dff217c3f237e37a *R/predict.radline.R
 194393ef5dc2d2ebb82d6641de089167 *R/predict.rda.R
 eb223fbfded71ae4f0b374c1e92c3f2e *R/predict.specaccum.R
 4f56d16f5bf8f9af3477c23137a70fb5 *R/pregraphKM.R
@@ -264,7 +268,6 @@ eb223fbfded71ae4f0b374c1e92c3f2e *R/predict.specaccum.R
 7db2fd99abc08bf5e1341e5b74fb4617 *R/prestonfit.R
 953d32321b6e12a30209a8cda78244c9 *R/print.CCorA.R
 3a1584c7d991683a61271fb2fc002b73 *R/print.MOStest.R
-a70e6320123f0206b106483e01fffac4 *R/print.adipart.R
 1e07dd6a9eefb1d0af31a4db98c94305 *R/print.adonis.R
 15b31674cb74df69467902853a9254d1 *R/print.anosim.R
 7c99a949e155f76c3a223741afb39836 *R/print.betadisper.R
@@ -275,7 +278,6 @@ f757010a0187dc81e9b844df25c58640 *R/print.cca.R
 65e888e34fa8a8e1d5b577fbadb3161a *R/print.envfit.R
 ff355b68b19f8d8c29917ca33d4e8b8d *R/print.factorfit.R
 53efa849e48c5c91a51f42282237253c *R/print.fisherfit.R
-aa617465dd91f58e0281c19e9ee5ce40 *R/print.hiersimu.R
 6da316510cb840a6a9dd8d31d0e205af *R/print.humpfit.R
 b31dbaa6493fdda1f865f95b3e889aab *R/print.isomap.R
 7965811709121f361245c75e4932b434 *R/print.mantel.R
@@ -284,13 +286,12 @@ f92fd82d10ce91e2cba2239e756e1332 *R/print.mantel.correlog.R
 7e9b861d872cad8928c5ccfa4de01fa3 *R/print.monoMDS.R
 8bd5bbb931a97ddada79e4552bd614b8 *R/print.mrpp.R
 946b3b708190211e9eb1acc94ffa102d *R/print.mso.R
-f7eb8693202b925972060373b33ef14c *R/print.multipart.R
 7c074bf7870cb4c306c6181769b28a19 *R/print.nestedchecker.R
 eed481e994c01ec4d7b443fb8cafad46 *R/print.nesteddisc.R
 91c6a9b43c8b045d11a4b8834d1c9d47 *R/print.nestedn0.R
 0f8e3935f95b67b96e066435743bbee0 *R/print.nestednodf.R
 2f1732fffc2fb487420a910a1d3f5971 *R/print.nestedtemp.R
-95506521c2fca3d03335f36b20287848 *R/print.oecosimu.R
+58402a74333f1bd193dc7706f930a1b9 *R/print.oecosimu.R
 faf2620b1fbaec410af7b6e3159510ce *R/print.permat.R
 575da3562c07c6324e84288ac603b011 *R/print.permutest.betadisper.R
 f0c12622e4f250aacca5b7fabe54cbd1 *R/print.permutest.cca.R
@@ -321,16 +322,16 @@ c80f3931f05ab3066dfe93b98e737856 *R/print.varpart234.R
 819af0297e5d0a907f7fa91319c67e96 *R/profile.MOStest.R
 8159854c33821cea4cb77e34d882d79e *R/profile.fisherfit.R
 2f6b69115ea549102dad9b1b22c88034 *R/profile.humpfit.R
-9946d595da0d63ba44a1e056d7dd764e *R/protest.R
+6b42822347758e8299fd4448a17b4a52 *R/protest.R
 9169bd797963b5b121684de528651170 *R/rad.lognormal.R
-c98f27d1d141b809df4e800328b21d51 *R/rad.null.R
-0931015d1df809bb64ef7ae192b493ce *R/rad.preempt.R
+b129148e6efbbe1c45482c93d66f959b *R/rad.null.R
+949aca6b5bb7954a91819b17e515e396 *R/rad.preempt.R
 5a7e143e1292681c3d9b1e7b1b792aa6 *R/rad.zipf.R
 6780818aadc7b8c92c8f9a26a84b7dc0 *R/rad.zipfbrot.R
-bb27d541dc69c5d50900acd4b945579c *R/radfit.R
-3914c0f96fbcbebfa4844ce5547f074a *R/radfit.data.frame.R
+23fd677c0c8d6f4a8e6c6444d2cc8601 *R/radfit.R
+235a5213c266f1b2101c767ad1528888 *R/radfit.data.frame.R
 2f6d8082f39540bbe7cf0e0cf3c666c9 *R/radfit.default.R
-d8081f1fe1a0b3fed75768d355d2c9eb *R/radlattice.R
+1241d5e79bd4e9a9fbeb6a4bf4c86023 *R/radlattice.R
 33e294ffd0c75b25d36d18d4a3ad9884 *R/rankindex.R
 0a7c015a3700b34f29cb6bb0b94f68d2 *R/rarecurve.R
 6be7a6edec618f464dd7c783eca371f0 *R/rarefy.R
@@ -441,52 +442,53 @@ c51905bd025ccea2737527b6fca4a081 *data/mite.pcnm.rda
 ee3c343418d7cf2e435028adf93205f1 *data/sipoo.rda
 f87df84297865b5faf31e232e97a0f94 *data/varechem.rda
 82153b3e47807b926b77cef49900f422 *data/varespec.rda
-fe89ca7a45cf17172c295e0637c45915 *inst/ChangeLog
-c1072e0d722947cb1f14848474f8a497 *inst/NEWS.Rd
+e277418485f322c350fb3f0f3060ca4c *inst/ChangeLog
+0929874442ebca450bb3dfb965e1a664 *inst/NEWS.Rd
 9abfab8b05c34dd283379a7d87500ffb *inst/ONEWS
-2bf5f9d0d6a982272e2f6a7a154b159b *inst/doc/FAQ-vegan.pdf
+276ec5d3b110d6c989a167df0be2ac60 *inst/doc/FAQ-vegan.pdf
 106b655f1410bb42eb5cb46b733451ec *inst/doc/FAQ-vegan.texi
-57c45e7d2bd1713cf12ab2b2bd20e0a4 *inst/doc/Makefile
-3c900692f99137362b10da3e15785d6b *inst/doc/NEWS.html
+a712f5f699b09276e523f83ce3b731ae *inst/doc/Makefile
+1f885635aeb01b38fae7c6a9178d0cf5 *inst/doc/NEWS.html
 bfad8283cb33e6e2837f62360671a014 *inst/doc/decision-vegan.Rnw
-6f974abe20a7b03315e5ff97db7bd6bf *inst/doc/decision-vegan.pdf
-96080f151551fe55e0fede9277229056 *inst/doc/decision-vegan.tex
-68b9e795bdeba8df6da627413bce041c *inst/doc/diversity-vegan.Rnw
-9b5980940f30caf9f03ee02f5d91d9b2 *inst/doc/diversity-vegan.pdf
-70b04c5745010d664b334bb08a74bdb5 *inst/doc/diversity-vegan.tex
-387e8ebdd98b55308b7d8e833a2c9c2a *inst/doc/intro-vegan.Rnw
-4da9c9eac1ef9715b588a3abfe8b93f4 *inst/doc/intro-vegan.pdf
-d402b6d714e9c382053ffae26ed8b5cb *inst/doc/intro-vegan.tex
+ed620c74dffdcb44f8e67777cc827e92 *inst/doc/decision-vegan.pdf
+95d0308a4a52046bd422a43d6fb89b10 *inst/doc/decision-vegan.tex
+5b8e069938a0131e3fef6efa8672c2f0 *inst/doc/diversity-vegan.Rnw
+bc98f37c4e876c1b40c44cb8dd6e09fa *inst/doc/diversity-vegan.pdf
+1a5efe2e0b2794c10cb6ae87cbde7e63 *inst/doc/diversity-vegan.tex
+81b0f1cea43478f0d5ed354b8eb8cb73 *inst/doc/intro-vegan.Rnw
+a0312443c08df5cc22f648b53bb20648 *inst/doc/intro-vegan.pdf
+8847f6fe515c51cfdf8c7cab1f14a77b *inst/doc/intro-vegan.tex
 a1c35ea488b715441cd2269eb6998945 *inst/doc/partitioning.pdf
 bafe698f5c543e20d4c52659545677bf *inst/doc/vegan.bib
 8790e1afcfbf88c374d6c87928eab54f *inst/doc/veganjss.sty
-aa69f7b97270f25d2a9871d33ff628bd *man/BCI.Rd
+e2d32a2a53e75e8be60574c7e1cc3239 *man/BCI.Rd
 33db614085aa448f4241cd79ddc62461 *man/CCorA.Rd
 525d1213753747626c919c22d14073f5 *man/MDSrotate.Rd
 fd218be03aa2591e5123d11780ccba1a *man/MOStest.Rd
-0c2b91d47a27d3eca28c1937172ad025 *man/RsquareAdj.Rd
+f2823a48acb6f861404b6682b3f52a45 *man/RsquareAdj.Rd
 eb6df0e2eba46d2bbdc056baff276481 *man/SSarrhenius.Rd
-52e5b417de274fcc3415c4a24ebfd167 *man/add1.cca.Rd
-40859386dd6bcf4f534086eb5af441be *man/adipart.Rd
-59a0c8a693e78a8bbfc697dc7e7c659e *man/adonis.Rd
-e53842f5c91b0d5a6405018455c054b9 *man/anosim.Rd
-4455b3a62d0afced71232987ce954345 *man/anova.cca.Rd
+d7499407f3dd1cc977ea094eaa2af35e *man/add1.cca.Rd
+f7a0445e698dc21782cd113d867f7142 *man/adipart.Rd
+88eaec1454840ff74dc06a22da6da672 *man/adonis.Rd
+0443d02a7f37af2e14e745c3c1eb3aee *man/anosim.Rd
+427787a2ec7a1cd49df31a08cfa1d74b *man/anova.cca.Rd
 c57af27fa11dadcd48981fcf42b2d221 *man/as.mlm.Rd
-c1fec65743ee33818457d6a8512f899b *man/beals.Rd
-b2067f585b93a5402e5e03e551f3974e *man/betadisper.Rd
-43d3d85bd4fb2968397cc8284dc28c9b *man/betadiver.Rd
-050375e9a158a5c1d01118dacb89553e *man/bgdispersal.Rd
-d70c076364daf2d719df7a8d30603d2d *man/bioenv.Rd
+8e3718248ff8d48e724654ab17caa2e2 *man/beals.Rd
+362efe66c82eee9e301f516cfb8ab25f *man/betadisper.Rd
+1336f0afb69a05bee9f6e7706d81d038 *man/betadiver.Rd
+b04c2fae35dba2d97cb248814d5e2fe9 *man/bgdispersal.Rd
+890f6ec54dcf7962cbaa197ef7be8339 *man/bioenv.Rd
 1eab4a6369fa1d203a4a3f41f4ee4c06 *man/biplot.rda.Rd
-b9c71f79ddf9ef942e4d66ea50077aae *man/capscale.Rd
-5bf9340fd1b1fe55841ceb14ceb47134 *man/cascadeKM.Rd
-4855513776fb88f9804527f2eacff70c *man/cca.Rd
-5a6355fecddeee09323325f7e8573f93 *man/cca.object.Rd
-0c376bfe24a62389460c7fa7d5d78b90 *man/clamtest.Rd
+e2cfaeba8e49d98db1952d2406d9b2ae *man/capscale.Rd
+d3c1067cb5e4dc6f6b9a1c7d94e70ab5 *man/cascadeKM.Rd
+3dd123c139b485271485355079c685a5 *man/cca.Rd
+46256ab1491bec3fe5cd98c0aaaebec0 *man/cca.object.Rd
+c083092c765a7249ba073c25c43c34a3 *man/clamtest.Rd
 335d0f7691ad9d0c48fffce9f9db6201 *man/contribdiv.Rd
 c41033fb9c572365490cc23b9870c950 *man/decorana.Rd
-47583461eb570b83ad1944b97668fad2 *man/decostand.Rd
-a62d730c922e7d476d45e8f21178f1b3 *man/designdist.Rd
+00a34002a7464d1008b2ba63526a4afe *man/decostand.Rd
+5a48fad137cfbbc09ff287edb5170542 *man/density.adonis.Rd
+22e3451a1cc9e294c2ad0e1a4531b136 *man/designdist.Rd
 c01e0664652fbc8ef4963059bee4e422 *man/deviance.cca.Rd
 fa5821f3d5077efdf01f78f68bcaec83 *man/dispindmorisita.Rd
 f3f742efa7511a4c33108a00b512ebd9 *man/distconnected.Rd
@@ -494,26 +496,26 @@ b9e1c696647b4029f944e830184eb2c9 *man/diversity.Rd
 31a227bd7e2dd4cf92c9b936e1f49963 *man/dune.Rd
 5bdeafda9c2d62feec3bde6b4cd49e3b *man/dune.taxon.Rd
 5f5f8c7df063606ccde6124c5dbe8add *man/eigenvals.Rd
-44da1163984ac12972ba0eb90aa50322 *man/envfit.Rd
+516c0b3d11d2746d15ead8919a35803c *man/envfit.Rd
 919976c144bda3a083af284255ea3695 *man/eventstar.Rd
-49360011c9efdee6d23c31bdc13e7061 *man/fisherfit.Rd
+5b7ba3db111582ad3e1f646bc8061ea8 *man/fisherfit.Rd
 841b3f32510ed2c3f64186d623f858ae *man/goodness.cca.Rd
-294c38ede8c6781d896d1beac104ad2b *man/goodness.metaMDS.Rd
-d1fa0f82d05b22f3607f372043022379 *man/humpfit.Rd
+4d5e44b51132481ab920292b2651041c *man/goodness.metaMDS.Rd
+81f199c3ba2c65a7b6f81cbb7cc9886d *man/humpfit.Rd
 c8fea575af3da292987d4f8c4aa831b0 *man/indpower.Rd
-3160a8b502bd11c8c195ec4df18159fd *man/isomap.Rd
-f05254c0a034b5131bc48c0642bec484 *man/kendall.global.Rd
+cebefbf4e3090e7f7605c2f5ce303da2 *man/isomap.Rd
+5f2c36639ab3e98a5fbb03bb8c8adda9 *man/kendall.global.Rd
 e473a6d2589993b85fc1176866fdde78 *man/linestack.Rd
-ecf0c9cbdf26bbadec0c055e618e94a6 *man/make.cepnames.Rd
-b2dea662ac33f79c7feed561c3e7ba5a *man/mantel.Rd
-9605f23ced71ce53c3a6a02c13e9a8af *man/mantel.correlog.Rd
-b8e7bb396124a40a0a9d6cca683139e2 *man/metaMDS.Rd
+a83cd418e2a052ea32ef431bce2e3787 *man/make.cepnames.Rd
+e8a27e8ace7bd56776794ebfd8391d64 *man/mantel.Rd
+f5db928e02639d251389edb1cd36ff13 *man/mantel.correlog.Rd
+8219156d6aa03c5882324feee963e5cb *man/metaMDS.Rd
 4cfb02239809fa03b28e10ec8e8c9c6b *man/mite.Rd
 0d0876437491b3701e4d20cbbf0d1b6f *man/model.matrix.cca.Rd
 54d0363b992caac0546642cd5c6d3f7b *man/monoMDS.Rd
-6c99ed635bd1c164710db16ded9e4402 *man/mrpp.Rd
+6ebc1e9be42fd2ffc21f4d8a81512f91 *man/mrpp.Rd
 99129fa5cefe7d1a83ca04f32f0a91ba *man/mso.Rd
-0b83da5291da2f3ee0fd387ef2bbca6f *man/multipart.Rd
+56a812e640082ccf928224114330130b *man/multipart.Rd
 17b0186a4e98a2fe6b66d9f56ad3f958 *man/nestedtemp.Rd
 2907c7465be1247d2ed930025b877403 *man/nobs.adonis.Rd
 ba28e7e71d0637fd962d53d1e23488c0 *man/oecosimu.Rd
@@ -521,11 +523,11 @@ ba28e7e71d0637fd962d53d1e23488c0 *man/oecosimu.Rd
 03aab4cb7ca71141281d2abd3e810231 *man/ordihull.Rd
 8f8a34c5fcfcc1fe9f88ca16e84a1da6 *man/ordilabel.Rd
 994cfc973f88c682b741e48377e1b9b4 *man/ordiplot.Rd
-10e145cf5d658345925da8536e2afb63 *man/ordiplot3d.Rd
+b88713c8f69a6a0e677d406b51d2062e *man/ordiplot3d.Rd
 ca2ed47e3adb6f4ab9aa861872b07e4b *man/ordipointlabel.Rd
 d4d27a34b2e7d9d1b732a0d06cb9d9f4 *man/ordiresids.Rd
 9bd611f72fad857e16a10f46f8f80e86 *man/ordistep.Rd
-1a3cb1337bb8ce725d62c4cb530a2f02 *man/ordisurf.Rd
+7b62153e759a450cfececacbd385cd63 *man/ordisurf.Rd
 42dfa56db8adbd0a44e716086dee6d7d *man/orditkplot.Rd
 f3d8ca5afe5c64965ac7cb658980d93f *man/orditorp.Rd
 d971701b3c6f89b3a6b358a3966a43d2 *man/ordixyplot.Rd
@@ -535,35 +537,35 @@ d971701b3c6f89b3a6b358a3966a43d2 *man/ordixyplot.Rd
 fb24c58ca61caf767ced5bd79dbff57e *man/permutest.betadisper.Rd
 a8862a670b6894fc38bbb2ade4c85e23 *man/plot.cca.Rd
 9f296f0830397598ee4849fc006762d1 *man/prc.Rd
-f72046fbb4bfbcbc13b60d4288f3e16c *man/predict.cca.Rd
-e0fcf2ebc4cfbad289fbc9d6dc4b74d2 *man/procrustes.Rd
+a02878ee16916bc9cdccb7bc6cd6b6c0 *man/predict.cca.Rd
+5208c960ebdfa5d51e6430763c8a944e *man/procrustes.Rd
 01a6ca946df5ad493adfb54003ad8e00 *man/pyrifos.Rd
-b8768707b4d121b6c4ea9feca9f472ce *man/radfit.Rd
+26f6e755b8b18407ed6c30dbf7887533 *man/radfit.Rd
 3e70bfa0a8ae5d4c3c60dba77500b584 *man/rankindex.Rd
-aa050f298d26952caa4343b8cfb3ef03 *man/raupcrick.Rd
+64342c9ea7e7b2607d433c3346f9726a *man/raupcrick.Rd
 2867f5f71a47da498cbadf9aaa01b2b6 *man/read.cep.Rd
 c1ffd9ca78ad968e1da11fdba007cbe8 *man/renyi.Rd
 5c25a88ca55fabce5783509c706faad5 *man/scores.Rd
-ce3aa1e8e15c5f95a345fc822630767b *man/screeplot.cca.Rd
+9732a76d9f971df9db16b97d5746615e *man/screeplot.cca.Rd
 947c357c856ef350340eb54673c0bc5c *man/simper.Rd
 45cd418b2264b4eb6abc89cc11a7877f *man/simulate.rda.Rd
 b34910fa6ed6c9bfbd90a7f7443a135f *man/sipoo.Rd
 d7dd63e022633049766cffdaf6cac723 *man/spantree.Rd
-48506cdf76f6d0ed87fd163c80e438eb *man/specaccum.Rd
+8d9d2fe132aaa595df4d7652dbe68c54 *man/specaccum.Rd
 53818a4edb1d52d425065bea76963021 *man/specpool.Rd
-756d6587c4bbcf81b568b0c02568f9e7 *man/stepacross.Rd
+5b9e51c85395f80f8504954e4175f877 *man/stepacross.Rd
 0aac5f5c8f58fc8fe1cb6c0ba819b196 *man/taxondive.Rd
 85f77fcf89b48586502c00baef8e5561 *man/tolerance.Rd
-11de9424624c14deb93f668009fd5992 *man/treedive.Rd
+af0e9dd8ea39743f9e7a8c471eed29e7 *man/treedive.Rd
 5a1c08a3f35d82027259b3000e94cd2e *man/tsallis.Rd
 033dd7d7917185cea81e4d7afcd59df9 *man/varechem.Rd
-e6e79b26209f4eeb4d8dab08a8698f6b *man/varpart.Rd
+be68944c43a32540e71a0e0b215be64b *man/varpart.Rd
 699122da39bdbbfbfeb6a1f8f078242c *man/vegan-defunct.Rd
 c8e9610be158d93af56d7db754e560d9 *man/vegan-deprecated.Rd
-7c69d0cfd6b295baf8eec54bbc99396c *man/vegan-internal.Rd
+e8d8b52de649979b17b484d03775c5c0 *man/vegan-internal.Rd
 ec9c7e972e815897f61fb3fdd16a008b *man/vegan-package.Rd
 f6f284ceb3b9a39e7b079936f9400cc2 *man/vegandocs.Rd
-72431d0c56876fa71afc972d9016d3fe *man/vegdist.Rd
+197fe8ed2fe66db294381cb3c08325e1 *man/vegdist.Rd
 3fc88b0b479d2099ced8ca764fbcbea6 *man/vegemite.Rd
 c3209a8eff0fe638d3a43b25ea5bec16 *man/wascores.Rd
 99d6eca6ce44322c03850991ef4de828 *man/wcmdscale.Rd
diff --git a/NAMESPACE b/NAMESPACE
index 7595e66..164286f 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -73,6 +73,7 @@ S3method(adipart, default)
 S3method(adipart, formula)
 # AIC: stats
 S3method(AIC, radfit)
+S3method(AIC, radfit.frame)
 # RsquareAdj: vegan
 S3method(RsquareAdj, cca)
 S3method(RsquareAdj, default)
@@ -123,6 +124,7 @@ S3method(cca, formula)
 # coef: stats
 S3method(coef, cca)
 S3method(coef, radfit)
+S3method(coef, radfit.frame)
 S3method(coef, rda)
 # confint: stats -- also uses MASS:::confint.glm & MASS:::profile.glm
 # does this work with namespaces??
@@ -131,12 +133,21 @@ S3method(confint, fisherfit)
 # cophenetic: stats
 S3method(cophenetic, spantree)
 # density: stats
+S3method(density, adonis)
+S3method(density, anosim)
+S3method(density, mantel)
+S3method(density, mrpp)
 S3method(density, oecosimu)
+S3method(density, permutest.cca)
+S3method(density, protest)
 # densityplot: lattice
+S3method(densityplot, adonis)
 S3method(densityplot, oecosimu)
 # deviance: stats
 S3method(deviance, cca)
 S3method(deviance, rda)
+S3method(deviance, radfit)
+S3method(deviance, radfit.frame)
 # drop1: stats
 S3method(drop1, cca)
 # eigenvals: vegan
@@ -163,6 +174,7 @@ S3method(fitted, capscale)
 S3method(fitted, cca)
 S3method(fitted, procrustes)
 S3method(fitted, radfit)
+S3method(fitted, radfit.frame)
 S3method(fitted, rda)
 # goodness: vegan
 S3method(goodness, cca)
@@ -176,6 +188,8 @@ S3method(hiersimu, default)
 S3method(hiersimu, formula)
 # identify: graphics
 S3method(identify, ordiplot)
+# labels: base
+S3method(labels, envfit)
 # lines: graphics
 S3method(lines, humpfit)
 S3method(lines, permat)
@@ -183,7 +197,11 @@ S3method(lines, preston)
 S3method(lines, prestonfit)
 S3method(lines, procrustes)
 S3method(lines, radline)
+S3method(lines, radfit)
 S3method(lines, spantree)
+## logLik: stats
+S3method(logLik, radfit)
+S3method(logLik, radfit.frame)
 # model.frame, model.matrix: stats
 S3method(model.frame, cca)
 S3method(model.matrix, cca)
@@ -254,6 +272,7 @@ S3method(plot, specaccum)
 S3method(plot, taxondive)
 S3method(plot, varpart)
 S3method(plot, varpart234)
+S3method(plot, vegandensity)
 # points: graphics
 S3method(points, cca)
 S3method(points, decorana)
@@ -263,18 +282,21 @@ S3method(points, ordiplot)
 S3method(points, orditkplot)
 S3method(points, procrustes)
 S3method(points, radline)
+S3method(points, radfit)
 # predict: stats
 S3method(predict, cca)
 S3method(predict, decorana)
 S3method(predict, fitspecaccum)
 S3method(predict, humpfit)
 S3method(predict, procrustes)
+S3method(predict, radline)
+S3method(predict, radfit)
+S3method(predict, radfit.frame)
 S3method(predict, rda)
 S3method(predict, specaccum)
 # print: base
 S3method(print, CCorA)
 S3method(print, MOStest)
-S3method(print, adipart)
 S3method(print, adonis)
 S3method(print, anosim)
 S3method(print, betadisper)
@@ -286,7 +308,6 @@ S3method(print, eigenvals)
 S3method(print, envfit)
 S3method(print, factorfit)
 S3method(print, fisherfit)
-S3method(print, hiersimu)
 S3method(print, humpfit)
 S3method(print, isomap)
 S3method(print, mantel)
@@ -295,7 +316,6 @@ S3method(print, metaMDS)
 S3method(print, monoMDS)
 S3method(print, mrpp)
 S3method(print, mso)
-S3method(print, multipart)
 S3method(print, nestedchecker)
 S3method(print, nesteddisc)
 S3method(print, nestedn0)
@@ -339,6 +359,7 @@ S3method(profile, humpfit)
 # radfit: vegan
 S3method(radfit, data.frame)
 S3method(radfit, default)
+S3method(radfit, matrix)
 # rda: vegan
 S3method(rda, default)
 S3method(rda, formula)
diff --git a/R/AIC.radfit.R b/R/AIC.radfit.R
index 04ac67a..dadd951 100644
--- a/R/AIC.radfit.R
+++ b/R/AIC.radfit.R
@@ -1,6 +1,41 @@
-"AIC.radfit" <-
+### these functions are defined _ex machina_ for radline objects which
+### inherit from glm. Here we define them for radfit objects where
+### object$models is a list of radline objects
+
+`AIC.radfit` <-
     function (object, k = 2, ...) 
 {
-    mods <- object$models
-    unlist(lapply(mods, AIC, k = k))
+    sapply(object$models, AIC, k = k, ...)
+}
+
+`deviance.radfit` <-
+    function(object, ...)
+{
+    sapply(object$models, deviance, ...)
+}
+
+`logLik.radfit` <-
+    function(object, ...)
+{
+    sapply(object$models, logLik, ...)
+}
+
+### Define also for radfit.frames which are lists of radfit objects
+
+`AIC.radfit.frame` <-
+    function(object, k = 2, ...)
+{
+    sapply(object, AIC, k = k, ...)
+}
+
+`deviance.radfit.frame` <-
+    function(object, ...)
+{
+    sapply(object, deviance, ...)
+}
+
+`logLik.radfit.frame` <-
+    function(object, ...)
+{
+    sapply(object, logLik, ...)
 }
diff --git a/R/adipart.default.R b/R/adipart.default.R
index 50c246d..d3b8c20 100644
--- a/R/adipart.default.R
+++ b/R/adipart.default.R
@@ -21,7 +21,7 @@ function(y, x, index=c("richness", "shannon", "simpson"),
         colnames(rhs) <- paste("level", 1:nlevs, sep="_")
     tlab <- colnames(rhs)
 
-    ## part check proper design of the model frame
+    ## check proper design of the model frame
     l1 <- sapply(rhs, function(z) length(unique(z)))
     if (!any(sapply(2:nlevs, function(z) l1[z] <= l1[z-1])))
         stop("number of levels are inapropriate, check sequence")
@@ -98,12 +98,14 @@ function(y, x, index=c("richness", "shannon", "simpson"),
     nam <- c(paste("alpha", 1:(nlevs-1), sep="."), "gamma",
              paste("beta", 1:(nlevs-1), sep="."))
     names(sim$statistic) <- attr(sim$oecosimu$statistic, "names") <- nam
-    attr(sim, "call") <- match.call()
-    attr(sim, "index") <- index
-    attr(sim, "weights") <- weights
+    call <- match.call()
+    call[[1]] <- as.name("adipart")
+    attr(sim, "call") <- call
+    attr(sim$oecosimu$simulated, "index") <- index
+    attr(sim$oecosimu$simulated, "weights") <- weights
     attr(sim, "n.levels") <- nlevs
     attr(sim, "terms") <- tlab
     attr(sim, "model") <- rhs
-    class(sim) <- c("adipart", "list")
+    class(sim) <- c("adipart", class(sim))
     sim
 }
diff --git a/R/adipart.formula.R b/R/adipart.formula.R
index 3d5c47c..e9f7265 100644
--- a/R/adipart.formula.R
+++ b/R/adipart.formula.R
@@ -1,113 +1,19 @@
-adipart.formula <-
-function(formula, data, index=c("richness", "shannon", "simpson"),
-    weights=c("unif", "prop"), relative = FALSE, nsimul=99, ...)
+`adipart.formula` <-
+    function(formula, data, index=c("richness", "shannon", "simpson"),
+             weights=c("unif", "prop"), relative = FALSE, nsimul=99, ...)
 {
     ## evaluate formula
-    lhs <- formula[[2]]
     if (missing(data))
         data <- parent.frame()
-    lhs <- as.matrix(eval(lhs, data))
-    formula[[2]] <- NULL
-    rhs <- model.frame(formula, data, drop.unused.levels = TRUE)
-    tlab <- attr(attr(rhs, "terms"), "term.labels")
-    nlevs <- length(tlab)
-    if (nlevs < 2)
-        stop("provide at least two level hierarchy")
-    if (any(rowSums(lhs) == 0))
-        stop("data matrix contains empty rows")
-    if (any(lhs < 0))
-        stop("data matrix contains negative entries")
+    tmp <- hierParseFormula(formula, data)
+    lhs <- tmp$lhs
+    rhs <- tmp$rhs
 
-    ## part check proper design of the model frame
-    noint <- attr(attr(attr(rhs, "terms"), "factors"), "dimnames")[[1]]
-    int <- attr(attr(attr(rhs, "terms"), "factors"), "dimnames")[[2]]
-    if (!identical(noint, int))
-        stop("interactions are not allowed in formula")
-    if (!all(attr(attr(rhs, "terms"), "dataClasses") == "factor"))
-        stop("all right hand side variables in formula must be factors")
-    l1 <- sapply(rhs, function(z) length(unique(z)))
-    if (!any(sapply(2:nlevs, function(z) l1[z] <= l1[z-1])))
-        stop("number of levels are inapropriate, check sequence")
-    rval <- list()
-    rval[[1]] <- as.factor(rhs[,nlevs])
-    rval[[1]] <- rval[[1]][drop = TRUE]
-    nCol <- nlevs - 1
-    for (i in 2:nlevs) {
-        rval[[i]] <- interaction(rhs[,nCol], rval[[(i-1)]], drop=TRUE)
-        nCol <- nCol - 1
-    }
-    rval <- as.data.frame(rval[rev(1:length(rval))])
-    l2 <- sapply(rval, function(z) length(unique(z)))
-    if (any(l1 != l2))
-        warning("levels are not perfectly nested")
-
-    ## aggregate response matrix
-    fullgamma <-if (nlevels(rhs[,nlevs]) == 1)
-        TRUE else FALSE
-    ftmp <- vector("list", nlevs)
-    for (i in 1:nlevs) {
-        ftmp[[i]] <- as.formula(paste("~", tlab[i], "- 1"))
-    }
-
-    ## is there a method/burnin/thin in ... ?
-    method <- if (is.null(list(...)$method))
-        "r2dtable" else list(...)$method
-    burnin <- if (is.null(list(...)$burnin))
-        0 else list(...)$burnin
-    thin <- if (is.null(list(...)$thin))
-        1 else list(...)$thin
-    base <- if (is.null(list(...)$base))
-        exp(1) else list(...)$base
-
-    ## evaluate other arguments
-    index <- match.arg(index)
-    weights <- match.arg(weights)
-    switch(index,
-           "richness" = {
-               divfun <- function(x) apply(x > 0, 1, sum)},
-           "shannon" = {
-               divfun <- function(x) diversity(x, index = "shannon", MARGIN = 1, base=base)},
-           "simpson" = {
-               divfun <- function(x) diversity(x, index = "simpson", MARGIN = 1)})
-
-    ## this is the function passed to oecosimu
-    wdivfun <- function(x) {
-        ## matrix sum *can* change in oecosimu (but default is constant sumMatr)
-        sumMatr <- sum(x)
-        if (fullgamma) {
-            tmp <- lapply(1:(nlevs-1), function(i) t(model.matrix(ftmp[[i]], rhs)) %*% x)
-            tmp[[nlevs]] <- matrix(colSums(x), nrow = 1, ncol = ncol(x))
-        } else {
-            tmp <- lapply(1:nlevs, function(i) t(model.matrix(ftmp[[i]], rhs)) %*% x)
-        }
-        ## weights will change in oecosimu thus need to be recalculated
-        if (weights == "prop")
-            wt <- lapply(1:nlevs, function(i) apply(tmp[[i]], 1, function(z) sum(z) / sumMatr))
-        else wt <- lapply(1:nlevs, function(i) rep(1 / NROW(tmp[[i]]), NROW(tmp[[i]])))
-        a <- sapply(1:nlevs, function(i) sum(divfun(tmp[[i]]) * wt[[i]]))
-        if (relative)
-            a <- a / a[length(a)]
-        b <- sapply(2:nlevs, function(i) a[i] - a[(i-1)])
-        c(a, b)
-    }
-    if (nsimul > 0) {
-        sim <- oecosimu(lhs, wdivfun, method = method, nsimul=nsimul,
-                        burnin=burnin, thin=thin)
-    } else {
-        sim <- wdivfun(lhs)
-        tmp <- rep(NA, length(sim))
-        sim <- list(statistic = sim,
-                    oecosimu = list(z = tmp, pval = tmp, method = NA, statistic = sim))
-    }
-    nam <- c(paste("alpha", 1:(nlevs-1), sep="."), "gamma",
-             paste("beta", 1:(nlevs-1), sep="."))
-    names(sim$statistic) <- attr(sim$oecosimu$statistic, "names") <- nam
-    attr(sim, "call") <- match.call()
-    attr(sim, "index") <- index
-    attr(sim, "weights") <- weights
-    attr(sim, "n.levels") <- nlevs
-    attr(sim, "terms") <- tlab
-    attr(sim, "model") <- rhs
-    class(sim) <- c("adipart", "list")
+    ## run simulations
+    sim <- adipart.default(lhs, rhs, index = index, weights = weights,
+                           relative = relative, nsimul = nsimul, ...)
+    call <- match.call()
+    call[[1]] <- as.name("adipart")
+    attr(sim, "call") <- call
     sim
 }
diff --git a/R/adonis.R b/R/adonis.R
index df19694..479d020 100644
--- a/R/adonis.R
+++ b/R/adonis.R
@@ -41,17 +41,18 @@
         dmat <- dist.lhs^2
     }
     n <- nrow(dmat)
-    I <- diag(n)
-    ones <- matrix(1,nrow=n)
-    A <- -(dmat)/2
-    G <- -.5 * dmat %*% (I - ones%*%t(ones)/n)
+    ## G is -dmat/2 centred by rows
+    G <- -sweep(dmat, 1, rowMeans(dmat))/2
     SS.Exp.comb <- sapply(H.s, function(hat) sum( G * t(hat)) )
     SS.Exp.each <- c(SS.Exp.comb - c(0,SS.Exp.comb[-nterms]) )
     H.snterm <- H.s[[nterms]]
+    ## t(I - H.snterm) is needed several times and we calculate it
+    ## here
+    tIH.snterm <- t(diag(n)-H.snterm)
     if (length(H.s) > 1)
         for (i in length(H.s):2)
             H.s[[i]] <- H.s[[i]] - H.s[[i-1]]
-    SS.Res <- sum( G * t(I-H.snterm))
+    SS.Res <- sum( G * tIH.snterm)
     df.Exp <- sapply(u.grps[-1], function(i) sum(grps==i) )
     df.Res <- n - qrhs$rank
     ## Get coefficients both for the species (if possible) and sites
@@ -66,8 +67,6 @@
     colnames(beta.sites) <- rownames(lhs)
     F.Mod <- (SS.Exp.each/df.Exp) / (SS.Res/df.Res)
 
-
-
     f.test <- function(tH, G, df.Exp, df.Res, tIH.snterm) {
       ## HERE I TRY CHANGING t(H)  TO tH, and
       ## t(I - H.snterm) to tIH.snterm, so that we don't have
@@ -98,7 +97,6 @@
 
 
     tH.s <- lapply(H.s, t)
-    tIH.snterm <- t(I-H.snterm)
     ## Apply permutations for each term
     ## This is the new f.test (2011-06-15) that uses fewer arguments
     f.perms <- sapply(1:nterms, function(i) {
diff --git a/R/anova.ccabyaxis.R b/R/anova.ccabyaxis.R
index 3e0748d..aff4b18 100644
--- a/R/anova.ccabyaxis.R
+++ b/R/anova.ccabyaxis.R
@@ -1,8 +1,6 @@
 `anova.ccabyaxis` <-
     function (object, cutoff = 1,  ...) 
 {
-    if(inherits(object, "prc"))
-        stop("anova(..., by = 'axis') cannot be used for 'prc' results")
     cutoff <- cutoff + sqrt(.Machine$double.eps)
     rnk <- object$CCA$rank
     if (!max(rnk, 0)) 
@@ -51,8 +49,10 @@
     environment(object$terms) <- environment()
     fla <- paste(". ~ ", axnam[1], "+ Condition(",
                  paste(axnam[-1], collapse="+"),")")
-    if (!is.null(CondMat))
+    if (!is.null(CondMat)) {
         fla <- paste(fla, " + Condition(CondMat)")
+        lc$CondMat <- CondMat
+    }
     fla <- update(formula(object), fla)
     sol <- anova(update(object, fla, data=lc),  ...)
     out[c(1, rnk + 1), ] <- sol
diff --git a/R/anova.ccabyterm.R b/R/anova.ccabyterm.R
index f60e987..e07a723 100644
--- a/R/anova.ccabyterm.R
+++ b/R/anova.ccabyterm.R
@@ -28,7 +28,7 @@
     if (!is.null(object$call$data))
         modelframe <- ordiGetData(object$call, globalenv())
     else
-        modelframe <- NULL
+        modelframe <- model.frame(object)
     for (.ITRM in ntrm:2) {
         if (ntrm < 2) 
             break
diff --git a/R/betadisper.R b/R/betadisper.R
index 9eb5a0a..1fc660c 100644
--- a/R/betadisper.R
+++ b/R/betadisper.R
@@ -1,6 +1,16 @@
 `betadisper` <-
     function(d, group, type = c("median","centroid"), bias.adjust=FALSE)
 {
+    ## inline function for double centring. We used .C("dblcen", ...,
+    ## PACKAGE = "stats") which does not dublicate its argument, but
+    ## it was removed from R in r60360 | ripley | 2012-08-22 07:59:00
+    ## UTC (Wed, 22 Aug 2012) "more conversion to .Call, clean up".
+    dblcen <- function(x, na.rm = TRUE) {
+        cnt <- colMeans(x, na.rm = na.rm)
+        x <- sweep(x, 2L, cnt, check.margin = FALSE)
+        cnt <- rowMeans(x, na.rm = na.rm)
+        sweep(x, 1L, cnt, check.margin = FALSE)
+    }
     ## inline function for spatial medians
     spatialMed <- function(vectors, group, pos) {
         axes <- seq_len(NCOL(vectors))
@@ -60,7 +70,7 @@
     }
     x <- x + t(x)
     storage.mode(x) <- "double"
-    .C("dblcen", x, as.integer(n), DUP = FALSE, PACKAGE="stats")
+    x <- dblcen(x)
     e <- eigen(-x/2, symmetric = TRUE)
     vectors <- e$vectors
     eig <- e$values
diff --git a/R/biplot.rda.R b/R/biplot.rda.R
index bb1d7ec..1f5fc11 100644
--- a/R/biplot.rda.R
+++ b/R/biplot.rda.R
@@ -23,18 +23,18 @@ biplot.rda <- function(x, choices = c(1, 2), scaling = 2,
   if (missing(type)) {
       nitlimit <- 80
       nit <- max(nrow(g$species), nrow(g$sites))
-      if (nit > nitlimit) 
-          type <- "points"
-      else type <- "text"
+      if (nit > nitlimit)
+          type <- rep("points", 2)
+      else type <- rep("text", 2)
   }
   else type <- match.arg(type, TYPES, several.ok = TRUE)
   if(length(type) < 2)
       type <- rep(type, 2)
-  if (missing(xlim)) 
+  if (missing(xlim))
       xlim <- range(g$species[, 1], g$sites[, 1])
-  if (missing(ylim)) 
+  if (missing(ylim))
       ylim <- range(g$species[, 2], g$sites[, 2])
-  plot(g[[1]], xlim = xlim, ylim = ylim, type = "n", asp = 1, 
+  plot(g[[1]], xlim = xlim, ylim = ylim, type = "n", asp = 1,
        ...)
   abline(h = 0, lty = 3)
   abline(v = 0, lty = 3)
@@ -51,9 +51,9 @@ biplot.rda <- function(x, choices = c(1, 2), scaling = 2,
                col = col[2], cex = 0.7)
   }
   if (!is.null(g$sites)) {
-      if (type[2] == "text") 
+      if (type[2] == "text")
           text(g$sites, rownames(g$sites), cex = 0.7, col = col[1])
-      else if (type[2] == "points") 
+      else if (type[2] == "points")
           points(g$sites, pch = 1, cex = 0.7, col = col[1])
   }
   class(g) <- "ordiplot"
diff --git a/R/cIndexKM.R b/R/cIndexKM.R
index 7301c44..d04bc8c 100644
--- a/R/cIndexKM.R
+++ b/R/cIndexKM.R
@@ -1,75 +1,7 @@
-"cIndexKM" <- function (y, x, index = "all") 
+`cIndexKM` <-
+    function (y, x, index = "all") 
 {
     kmeans_res <- y
-#########################################
-    withinss <- function(kmeans_res, x) 
-    {
-        retval <- rep(0, nrow(kmeans_res$centers))
-        x <- (x - kmeans_res$centers[kmeans_res$cluster, ])^2
-        for (k in 1:nrow(kmeans_res$centers)) 
-        {
-            retval[k] <- sum(x[kmeans_res$cluster == k, ])
-        }
-        retval
-    }
-##########################################
-    varwithinss <- function(x, centers, cluster) 
-    {
-        nrow <- dim(centers)[1]
-        nvar <- dim(x)[2]
-        varwithins <- matrix(0, nrow, nvar)
-        x <- (x - centers[cluster, ])^2
-        for (l in 1:nvar) 
-        {
-            for (k in 1:nrow) 
-            {
-                varwithins[k, l] <- sum(x[cluster == k, l])
-            }
-        }
-        varwithins
-    }
-##########################################
-    maxmindist <- function(clsize, distscen) 
-    {
-        ncl <- length(clsize)
-        npairs <- 0
-        for (i in 1:ncl) npairs <- npairs + clsize[i] * (clsize[i] - 1)/2
-        mindw <- 0
-        nfound <- distscen[1]
-        i <- 1
-        while (nfound < npairs) 
-        {
-            if ((nfound + distscen[i + 1]) < npairs) 
-            {
-                mindw <- mindw + i * distscen[i + 1]
-                nfound <- nfound + distscen[i + 1]
-            }
-            else 
-            {
-                mindw <- mindw + i * (npairs - nfound)
-                nfound <- nfound + distscen[i + 1]
-            }
-            i <- i + 1
-        }
-        maxdw <- 0
-        nfound <- 0
-        i <- length(distscen) - 1
-        while (nfound < npairs) 
-        {
-            if ((nfound + distscen[i + 1]) < npairs) 
-            {
-                maxdw <- maxdw + i * distscen[i + 1]
-                nfound <- nfound + distscen[i + 1]
-            }
-            else 
-            {
-                maxdw <- maxdw + i * (npairs - nfound)
-                nfound <- nfound + distscen[i + 1]
-            }
-            i <- i - 1
-        }
-        list(mindw = mindw, maxdw = maxdw)
-    }
 #############################################
     gss <- function(x, clsize, withins) 
     {
@@ -83,39 +15,7 @@
         list(wgss = wgss, bgss = bgss)
     }
 #############################################
-    vargss <- function(x, clsize, varwithins) 
-    {
-        nvar <- dim(x)[2]
-        n <- sum(clsize)
-        k <- length(clsize)
-        varallmean <- rep(0, nvar)
-        varallmeandist <- rep(0, nvar)
-        varwgss <- rep(0, nvar)
-        for (l in 1:nvar) varallmean[l] <- mean(x[, l])
-        vardmean <- sweep(x, 2, varallmean, "-")
-        for (l in 1:nvar) 
-        {
-            varallmeandist[l] <- sum((vardmean[, l])^2)
-            varwgss[l] <- sum(varwithins[, l])
-        }
-        varbgss <- varallmeandist - varwgss
-        vartss <- varbgss + varwgss
-        list(vartss = vartss, varbgss = varbgss)
-    }
-		
-#################################################
-    count <- function(x) 
-    {
-        nr <- nrow(x)
-        nc <- ncol(x)
-        d <- integer(nc + 1)
-        retval <- .C("count", xrows = nr, xcols = nc, x = as.integer(x), 
-                     d = d, PACKAGE = "cclust")
-        d <- retval$d
-        names(d) <- 0:nc
-        d
-    }
-################################################
+
 ### Function modified by SD and PL from the original "cIndexKM" in "cclust" 
 ### to accommodate a single response variable as well as singleton groups
 ### and remove unwanted index.
diff --git a/R/coef.radfit.R b/R/coef.radfit.R
index 2398246..9e121fd 100644
--- a/R/coef.radfit.R
+++ b/R/coef.radfit.R
@@ -1,4 +1,4 @@
-"coef.radfit" <-
+`coef.radfit` <-
     function (object, ...) 
 {
     out <- sapply(object$models, function(x) if (length(coef(x)) < 
@@ -9,3 +9,9 @@
     colnames(out) <- paste("par", 1:3, sep = "")
     out
 }
+
+`coef.radfit.frame` <-
+    function(object, ...)
+{
+    lapply(object, coef, ...)
+}
diff --git a/R/density.anosim.R b/R/density.anosim.R
new file mode 100644
index 0000000..4dff5bb
--- /dev/null
+++ b/R/density.anosim.R
@@ -0,0 +1,136 @@
+### density & densityplot methods for vegan functions returning
+### statistics from permuted/simulated data. These are modelled after
+### density.oecosimu and densityplot.oecosimu (which are in their
+### separate files).
+
+## anosim
+
+`density.anosim` <-
+    function(x, ...)
+{
+    obs <- x$statistic
+    ## Put observed statistic among permutations
+    out <- density(c(obs, x$perm), ...)
+    out$call <- match.call()
+    out$observed <- obs
+    out$call[[1]] <- as.name("density")
+    class(out) <- c("vegandensity", class(out))
+    out
+}
+
+## adonis can return a matrix of terms, hence we also have densityplot()
+
+`density.adonis` <-
+    function(x, ...)
+{
+    cols <- ncol(x$f.perms)
+    if (cols > 1)
+        warning("'density' is meaningful only with one term, you have ", cols)
+    obs <- x$aov.tab$F.Model
+    obs <- obs[!is.na(obs)]
+    out <- density(c(obs, x$f.perms), ...)
+    out$observed <- obs
+    out$call <- match.call()
+    out$call[[1]] <- as.name("density")
+    class(out) <- c("vegandensity", class(out))
+    out
+}
+
+`densityplot.adonis` <-
+    function(x, data, xlab = "Null", ...)
+{
+    require(lattice) || stop("requires package 'lattice'")
+    obs <- x$aov.tab$F.Model
+    obs <- obs[!is.na(obs)]
+    sim <- rbind(obs, x$f.perms)
+    nm <- rownames(x$aov.tab)[col(sim)]
+    densityplot( ~ as.vector(sim) | factor(nm, levels = unique(nm)),
+                xlab = xlab,
+                panel = function(x, ...) {
+                    panel.densityplot(x, ...)
+                    panel.abline(v = obs[panel.number()], ...)
+                },
+                ...)
+}
+
+## mantel
+
+`density.mantel` <-
+    function(x, ...)
+{
+    obs <- x$statistic
+    out <- density(c(obs, x$perm), ...)
+    out$observed <- obs
+    out$call <- match.call()
+    out$call[[1]] <- as.name("density")
+    class(out) <- c("vegandensity", class(out))
+    out
+}
+
+## mrpp
+
+`density.mrpp` <-
+    function(x, ...)
+{
+    obs <- x$delta
+    out <- density(c(obs, x$boot.deltas), ...)
+    out$observed <- obs
+    out$call <- match.call()
+    out$call[[1]] <- as.name("density")
+    class(out) <- c("vegandensity", class(out))
+    out
+}
+
+## anova.cca does not return permutation results, but permutest.cca
+## does. However, permutest.cca always finds only one statistic. Full
+## tables anova.cca are found by repeated calls to permutest.cca.
+
+`density.permutest.cca` <-
+    function(x, ...)
+{
+    obs <- x$F.0
+    out <- density(c(obs, x$F.perm), ...)
+    out$observed <- obs
+    out$call <- match.call()
+    out$call[[1]] <- as.name("density")
+    class(out) <- c("vegandensity", class(out))
+    out
+}
+
+## protest
+
+`density.protest` <-
+    function(x, ...)
+{
+    obs <- x$t0
+    out <- density(c(obs, x$t), ...)
+    out$observed <- obs
+    out$call <- match.call()
+    out$call[[1]] <- as.name("density")
+    class(out) <- c("vegandensity", class(out))
+    out
+}
+
+#### plot method: the following copies stats::plot.density() code but
+#### adds one new argument to draw abline(v=...) for the observed
+#### statistic
+
+`plot.vegandensity` <-
+    function (x, main = NULL, xlab = NULL, ylab = "Density", type = "l", 
+    zero.line = TRUE, obs.line = TRUE, ...) 
+{
+    if (is.null(xlab)) 
+        xlab <- paste("N =", x$n, "  Bandwidth =", formatC(x$bw))
+    if (is.null(main)) 
+        main <- deparse(x$call)
+    ## change obs.line to col=2 (red) if it was logical TRUE
+    if (isTRUE(obs.line))
+        obs.line <- 2
+    plot.default(x, main = main, xlab = xlab, ylab = ylab, type = type,
+                 ...)
+    if (zero.line) 
+        abline(h = 0, lwd = 0.1, col = "gray")
+    if (is.character(obs.line) || obs.line)
+        abline(v = x$observed, col = obs.line)
+    invisible(NULL)
+}
diff --git a/R/density.oecosimu.R b/R/density.oecosimu.R
index 25e01ca..640df4c 100644
--- a/R/density.oecosimu.R
+++ b/R/density.oecosimu.R
@@ -1,8 +1,14 @@
 `density.oecosimu` <-
     function(x, ...)
 {
-    out <- density(t(x$oecosimu$simulated), ...)
+    cols <- nrow(x$oecosimu$simulated)
+    if (cols > 1)
+        warning("'density' is meaningful only with one statistic, you have ", cols)
+    obs <- x$oecosimu$statistic
+    out <- density(rbind(obs, t(x$oecosimu$simulated)), ...)
+    out$observed <- obs
     out$call <- match.call()
+    out$call[[1]] <- as.name("density")
+    class(out) <- c("vegandensity", class(out))
     out
 }
-
diff --git a/R/densityplot.oecosimu.R b/R/densityplot.oecosimu.R
index ba609b6..2700fe5 100644
--- a/R/densityplot.oecosimu.R
+++ b/R/densityplot.oecosimu.R
@@ -2,8 +2,8 @@
     function(x, data, xlab = "Simulated", ...)
 {
     require(lattice) || stop("requires package 'lattice'")
-    sim <- t(x$oecosimu$simulated)
     obs <- x$oecosimu$statistic
+    sim <- rbind(obs, t(x$oecosimu$simulated))
     nm <- names(obs)[col(sim)]
     densityplot( ~ as.vector(sim) | factor(nm, levels = unique(nm)),
                 xlab = xlab,
diff --git a/R/fitted.radfit.R b/R/fitted.radfit.R
index 3db59d5..958f4a7 100644
--- a/R/fitted.radfit.R
+++ b/R/fitted.radfit.R
@@ -1,5 +1,11 @@
 `fitted.radfit` <-
     function(object, ...)
 {
-    matrix(sapply(object$models, fitted), ncol=length(object$models))
+    sapply(object$models, fitted)
+}
+
+`fitted.radfit.frame` <-
+    function(object, ...)
+{
+    lapply(object, fitted, ...)
 }
diff --git a/R/hierParseFormula.R b/R/hierParseFormula.R
new file mode 100644
index 0000000..0d8e911
--- /dev/null
+++ b/R/hierParseFormula.R
@@ -0,0 +1,20 @@
+"hierParseFormula" <-
+function (formula, data)
+{
+    lhs <- formula[[2]]
+    if (any(attr(terms(formula, data = data), "order") > 1)) 
+        stop("interactions are not allowed")
+    lhs <- as.matrix(eval(lhs, data))
+    formula[[2]] <- NULL
+    rhs <- model.frame(formula, data, drop.unused.levels = TRUE)
+    rhs[] <- lapply(rhs, function(u) {
+        if (!is.factor(u)) 
+            u <- factor(u)
+        u
+    })
+    if (length(rhs) < 2)
+        stop("at least 2 hierarchy levels are needed")
+    attr(rhs, "terms") <- NULL
+    list(lhs=lhs, rhs=rhs)
+}
+
diff --git a/R/hiersimu.default.R b/R/hiersimu.default.R
index a3b29a0..621bbc0 100644
--- a/R/hiersimu.default.R
+++ b/R/hiersimu.default.R
@@ -15,7 +15,7 @@ relative = FALSE, drop.highest = FALSE, nsimul=99, ...)
         colnames(rhs) <- paste("level", 1:nlevs, sep="_")
     tlab <- colnames(rhs)
 
-    ## part check proper design of the model frame
+    ## check proper design of the model frame
     l1 <- sapply(rhs, function(z) length(unique(z)))
     if (!any(sapply(2:nlevs, function(z) l1[z] <= l1[z-1])))
         stop("number of levels are inapropriate, check sequence")
@@ -82,12 +82,14 @@ relative = FALSE, drop.highest = FALSE, nsimul=99, ...)
 #    nam <- paste("level", 1:nlevs, sep=".")
 #    names(sim$statistic) <- attr(sim$oecosimu$statistic, "names") <- nam
     names(sim$statistic) <- attr(sim$oecosimu$statistic, "names") <- tlab[1:nlevs]
-    attr(sim, "call") <- match.call()
+    call <- match.call()
+    call[[1]] <- as.name("hiersimu")
+    attr(sim, "call") <- call
     attr(sim, "FUN") <- FUN
     attr(sim, "location") <- location
     attr(sim, "n.levels") <- nlevs
     attr(sim, "terms") <- tlab
     attr(sim, "model") <- rhs
-    class(sim) <- c("hiersimu", "list")
+    class(sim) <- c("hiersimu", class(sim))
     sim
 }
diff --git a/R/hiersimu.formula.R b/R/hiersimu.formula.R
index 01eaaae..a2b9b11 100644
--- a/R/hiersimu.formula.R
+++ b/R/hiersimu.formula.R
@@ -1,96 +1,21 @@
-hiersimu.formula <-
-function(formula, data, FUN, location = c("mean", "median"),
-relative = FALSE, drop.highest = FALSE, nsimul=99, ...)
+`hiersimu.formula` <-
+    function(formula, data, FUN, location = c("mean", "median"),
+             relative = FALSE, drop.highest = FALSE, nsimul=99, ...)
 {
     ## evaluate formula
-    lhs <- formula[[2]]
     if (missing(data))
         data <- parent.frame()
-    lhs <- as.matrix(eval(lhs, data))
-    formula[[2]] <- NULL
-    rhs <- model.frame(formula, data, drop.unused.levels = TRUE)
-    tlab <- attr(attr(rhs, "terms"), "term.labels")
-    nlevs <- length(tlab)
+    tmp <- hierParseFormula(formula, data)
+    lhs <- tmp$lhs
+    rhs <- tmp$rhs
 
-    ## part check proper design of the model frame
-    noint <- attr(attr(attr(rhs, "terms"), "factors"), "dimnames")[[1]]
-    int <- attr(attr(attr(rhs, "terms"), "factors"), "dimnames")[[2]]
-    if (!identical(noint, int))
-        stop("interactions are not allowed in formula")
-    if (!all(attr(attr(rhs, "terms"), "dataClasses") == "factor"))
-        stop("all right hand side variables in formula must be factors")
-    l1 <- sapply(rhs, function(z) length(unique(z)))
-    if (nlevs > 1 && !any(sapply(2:nlevs, function(z) l1[z] <= l1[z-1])))
-        stop("number of levels are inapropriate, check sequence")
-    rval <- list()
-    rval[[1]] <- as.factor(rhs[,nlevs])
-    rval[[1]] <- rval[[1]][drop = TRUE]
-    if (nlevs > 1) {
-        nCol <- nlevs - 1
-        for (i in 2:nlevs) {
-            rval[[i]] <- interaction(rhs[,nCol], rval[[(i-1)]], drop=TRUE)
-            nCol <- nCol - 1
-        }
-    }
-    rval <- as.data.frame(rval[rev(1:length(rval))])
-    l2 <- sapply(rval, function(z) length(unique(z)))
-    if (any(l1 != l2))
-        warning("levels are not perfectly nested")
-
-    ## aggregate response matrix
-    fullgamma <-if (nlevels(rhs[,nlevs]) == 1)
-        TRUE else FALSE
-    if (fullgamma && drop.highest)
-        nlevs <- nlevs - 1
-    if (nlevs == 1 && relative)
-        stop("'relative=FALSE' makes no sense with 1 level")
-    ftmp <- vector("list", nlevs)
-    for (i in 1:nlevs) {
-        ftmp[[i]] <- as.formula(paste("~", tlab[i], "- 1"))
-    }
-
-    ## is there a method/burnin/thin in ... ?
-    method <- if (is.null(list(...)$method))
-        "r2dtable" else list(...)$method
-    burnin <- if (is.null(list(...)$burnin))
-        0 else list(...)$burnin
-    thin <- if (is.null(list(...)$thin))
-        1 else list(...)$thin
-
-    ## evaluate other arguments
-    if (!is.function(FUN))
-        stop("'FUN' must be a function")
-    location <- match.arg(location)
-    aggrFUN <- switch(location,
-        "mean" = mean,
-        "median" = median)
-
-    ## this is the function passed to oecosimu
-    evalFUN <- function(x) {
-        if (fullgamma && !drop.highest) {
-            tmp <- lapply(1:(nlevs-1), function(i) t(model.matrix(ftmp[[i]], rhs)) %*% x)
-            tmp[[nlevs]] <- matrix(colSums(x), nrow = 1, ncol = ncol(x))
-        } else {
-            tmp <- lapply(1:nlevs, function(i) t(model.matrix(ftmp[[i]], rhs)) %*% x)
-        }
-        a <- sapply(1:nlevs, function(i) aggrFUN(FUN(tmp[[i]]))) # dots removed from FUN
-        if (relative)
-            a <- a / a[length(a)]
-        a
-    }
-
-    ## processing oecosimu results
-    sim <- oecosimu(lhs, evalFUN, method = method, nsimul=nsimul,
-        burnin=burnin, thin=thin)
-#    nam <- paste("level", 1:nlevs, sep=".")
-#    names(sim$statistic) <- attr(sim$oecosimu$statistic, "names") <- nam
-    names(sim$statistic) <- attr(sim$oecosimu$statistic, "names") <- tlab[1:nlevs]
-    attr(sim, "call") <- match.call()
-    attr(sim, "FUN") <- FUN
-    attr(sim, "location") <- location
-    attr(sim, "n.levels") <- nlevs
-    attr(sim, "terms") <- tlab
-    attr(sim, "model") <- rhs
-    class(sim) <- c("hiersimu", "list")
+    ## run simulations
+    sim <- hiersimu.default(lhs, rhs, FUN = FUN, location = location,
+                            relative = relative, drop.highest = drop.highest,
+                            nsimul = nsimul, ...)
+    call <- match.call()
+    call[[1]] <- as.name("hiersimu")
+    attr(sim, "call") <- call
     sim
 }
+
diff --git a/R/labels.envfit.R b/R/labels.envfit.R
new file mode 100644
index 0000000..fd402d7
--- /dev/null
+++ b/R/labels.envfit.R
@@ -0,0 +1,9 @@
+`labels.envfit` <-
+    function(object, ...)
+{
+    out <- list("vectors" = rownames(object$vectors$arrows),
+                "factors" = rownames(object$factors$centroids))
+    if (is.null(out$vectors) || is.null(out$factors))
+        out <- unlist(out, use.names = FALSE)
+    out
+}
diff --git a/R/lines.radline.R b/R/lines.radline.R
index eda365e..85e0c4d 100644
--- a/R/lines.radline.R
+++ b/R/lines.radline.R
@@ -1,4 +1,4 @@
-"lines.radline" <-
+`lines.radline` <-
     function (x, ...) 
 {
     lin <- fitted(x)
@@ -6,3 +6,9 @@
     lines(rnk, lin, ...)
     invisible()
 }
+
+`lines.radfit` <-
+    function(x, ...)
+{
+    matlines(fitted(x), ...)
+}
diff --git a/R/mantel.R b/R/mantel.R
index 5d672bb..f0de5d8 100644
--- a/R/mantel.R
+++ b/R/mantel.R
@@ -1,24 +1,46 @@
 "mantel" <-
   function (xdis, ydis, method = "pearson", permutations = 999, 
-            strata) 
+            strata, na.rm = FALSE) 
 {
     xdis <- as.dist(xdis)
     ydis <- as.vector(as.dist(ydis))
-    tmp <- cor.test(as.vector(xdis), ydis, method = method)
-    statistic <- as.numeric(tmp$estimate)
-    variant <- tmp$method
+    ## Handle missing values
+    if (na.rm)
+        use <- "complete.obs"
+    else
+        use <- "all.obs"
+    statistic <- cor(as.vector(xdis), ydis, method = method, use = use)
+    variant <- match.arg(method, eval(formals(cor)$method))
+    variant <- switch(variant,
+                      pearson = "Pearson's product-moment correlation",
+                      kendall = "Kendall's rank correlation tau",
+                      spearman = "Spearman's rank correlation rho",
+                      variant)
+    N <- attr(xdis, "Size")
+    if (length(permutations) == 1) {
+        if (permutations > 0) {
+            arg <- if (missing(strata)) NULL else strata
+            permat <- t(replicate(permutations,
+                                  permuted.index(N, strata = arg)))
+        }
+    } else {
+        permat <- as.matrix(permutations)
+        if (ncol(permat) != N)
+            stop(gettextf("'permutations' have %d columns, but data have %d observations",
+                          ncol(permat), N))
+        permutations <- nrow(permutations)
+    }
     if (permutations) {
-        N <- attr(xdis, "Size")
-        perm <- rep(0, permutations)
-        ## asdist asn an index selects lower diagonal like as.dist,
-        ## but is faster since it does not seet 'dist' attributes
+        perm <- numeric(permutations)
+        ## asdist as an index selects lower diagonal like as.dist,
+        ## but is faster since it does not set 'dist' attributes
         xmat <- as.matrix(xdis)
         asdist <- row(xmat) > col(xmat)
-        for (i in 1:permutations) {
-            take <- permuted.index(N, strata)
+        ptest <- function(take, ...) {
             permvec <- (xmat[take, take])[asdist]
-            perm[i] <- cor(permvec, ydis, method = method)
+            drop(cor(permvec, ydis, method = method, use = use))
         }
+        perm <- sapply(1:permutations, function(i, ...) ptest(permat[i,], ...) )
         signif <- (sum(perm >= statistic) + 1)/(permutations + 1)
      }
     else {
diff --git a/R/mantel.partial.R b/R/mantel.partial.R
index 400e42a..eee21e5 100644
--- a/R/mantel.partial.R
+++ b/R/mantel.partial.R
@@ -1,6 +1,6 @@
 "mantel.partial" <-
   function (xdis, ydis, zdis, method = "pearson", permutations = 999, 
-            strata) 
+            strata, na.rm = FALSE) 
 {
     part.cor <- function(rxy, rxz, ryz) {
         (rxy - rxz * ryz)/sqrt(1-rxz*rxz)/sqrt(1-ryz*ryz)
@@ -8,11 +8,20 @@
     xdis <- as.dist(xdis)
     ydis <- as.vector(as.dist(ydis))
     zdis <- as.vector(as.dist(zdis))
-    rxy <- cor.test(as.vector(xdis), ydis, method = method)
-    rxz <- cor(as.vector(xdis), zdis, method = method)
-    ryz <- cor(ydis, zdis, method = method)
-    variant <- rxy$method
-    rxy <- rxy$estimate
+    ## Handle missing values
+    if (na.rm)
+        use <- "complete.obs"
+    else
+        use <- "all.obs"
+    rxy <- cor(as.vector(xdis), ydis, method = method, use = use)
+    rxz <- cor(as.vector(xdis), zdis, method = method, use = use)
+    ryz <- cor(ydis, zdis, method = method, use = use)
+    variant <- match.arg(method, eval(formals(cor)$method))
+    variant <- switch(variant,
+                      pearson = "Pearson's product-moment correlation",
+                      kendall = "Kendall's rank correlation tau",
+                      spearman = "Spearman's rank correlation rho",
+                      variant)
     statistic <- part.cor(rxy, rxz, ryz)
     if (permutations) {
         N <- attr(xdis, "Size")
@@ -22,8 +31,8 @@
         for (i in 1:permutations) {
             take <- permuted.index(N, strata)
             permvec <- (xmat[take, take])[asdist]
-            rxy <- cor(permvec, ydis, method = method)
-            rxz <- cor(permvec, zdis, method = method)
+            rxy <- cor(permvec, ydis, method = method, use = use)
+            rxz <- cor(permvec, zdis, method = method, use = use)
             perm[i] <- part.cor(rxy, rxz, ryz)
         }
         signif <- (sum(perm >= statistic)+1)/(permutations + 1)
diff --git a/R/multipart.default.R b/R/multipart.default.R
index c9d2c83..57fe78f 100644
--- a/R/multipart.default.R
+++ b/R/multipart.default.R
@@ -1,6 +1,6 @@
-multipart.default <-
-function(y, x, index=c("renyi", "tsallis"), scales = 1,
-    global = FALSE, relative = FALSE, nsimul=99, ...)
+`multipart.default` <-
+    function(y, x, index=c("renyi", "tsallis"), scales = 1,
+             global = FALSE, relative = FALSE, nsimul=99, ...)
 {
     if (length(scales) > 1)
         stop("length of 'scales' must be 1")
@@ -23,7 +23,7 @@ function(y, x, index=c("renyi", "tsallis"), scales = 1,
         colnames(rhs) <- paste("level", 1:nlevs, sep="_")
     tlab <- colnames(rhs)
 
-     ## part check proper design of the model frame
+     ## check proper design of the model frame
     l1 <- sapply(rhs, function(z) length(unique(z)))
     if (!any(sapply(2:nlevs, function(z) l1[z] <= l1[z-1])))
         stop("number of levels are inapropriate, check sequence")
@@ -123,13 +123,15 @@ function(y, x, index=c("renyi", "tsallis"), scales = 1,
     nam <- c(paste("alpha", 1:(nlevs-1), sep="."), "gamma",
         paste("beta", 1:(nlevs-1), sep="."))
     names(sim$statistic) <- attr(sim$oecosimu$statistic, "names") <- nam
-    attr(sim, "call") <- match.call()
-    attr(sim, "index") <- index
-    attr(sim, "scales") <- scales
-    attr(sim, "global") <- TRUE
+    call <- match.call()
+    call[[1]] <- as.name("multipart")
+    attr(sim, "call") <- call
+    attr(sim$oecosimu$simulated, "index") <- index
+    attr(sim$oecosimu$simulated, "scales") <- scales
+    attr(sim$oecosimu$simulated, "global") <- TRUE
     attr(sim, "n.levels") <- nlevs
     attr(sim, "terms") <- tlab
     attr(sim, "model") <- rhs
-    class(sim) <- c("multipart", "list")
+    class(sim) <- c("multipart", class(sim))
     sim
 }
diff --git a/R/multipart.formula.R b/R/multipart.formula.R
index b901aa8..e6c089e 100644
--- a/R/multipart.formula.R
+++ b/R/multipart.formula.R
@@ -1,139 +1,20 @@
-multipart.formula <-
-function(formula, data, index=c("renyi", "tsallis"), scales = 1,
-    global = FALSE, relative = FALSE, nsimul=99, ...)
+`multipart.formula` <-
+    function(formula, data, index=c("renyi", "tsallis"), scales = 1,
+             global = FALSE, relative = FALSE, nsimul=99, ...)
 {
-    if (length(scales) > 1)
-        stop("length of 'scales' must be 1")
     ## evaluate formula
-    lhs <- formula[[2]]
     if (missing(data))
         data <- parent.frame()
-    lhs <- as.matrix(eval(lhs, data))
-    formula[[2]] <- NULL
-    rhs <- model.frame(formula, data, drop.unused.levels = TRUE)
-    tlab <- attr(attr(rhs, "terms"), "term.labels")
-    nlevs <- length(tlab)
-    if (nlevs < 2)
-        stop("provide at least two level hierarchy")
-    if (any(rowSums(lhs) == 0))
-        stop("data matrix contains empty rows")
-    if (any(lhs < 0))
-        stop("data matrix contains negative entries")
+    tmp <- hierParseFormula(formula, data)
+    lhs <- tmp$lhs
+    rhs <- tmp$rhs
 
-    ## part check proper design of the model frame
-    noint <- attr(attr(attr(rhs, "terms"), "factors"), "dimnames")[[1]]
-    int <- attr(attr(attr(rhs, "terms"), "factors"), "dimnames")[[2]]
-    if (!identical(noint, int))
-        stop("interactions are not allowed in formula")
-    if (!all(attr(attr(rhs, "terms"), "dataClasses") == "factor"))
-        stop("all right hand side variables in formula must be factors")
-    l1 <- sapply(rhs, function(z) length(unique(z)))
-    if (!any(sapply(2:nlevs, function(z) l1[z] <= l1[z-1])))
-        stop("number of levels are inapropriate, check sequence")
-    rval <- list()
-    rval[[1]] <- as.factor(rhs[,nlevs])
-    rval[[1]] <- rval[[1]][drop = TRUE]
-    nCol <- nlevs - 1
-    for (i in 2:nlevs) {
-        rval[[i]] <- interaction(rhs[,nCol], rval[[(i-1)]], drop=TRUE)
-        nCol <- nCol - 1
-    }
-    rval <- as.data.frame(rval[rev(1:length(rval))])
-    l2 <- sapply(rval, function(z) length(unique(z)))
-    if (any(l1 != l2))
-        warning("levels are not perfectly nested")
-
-    ## aggregate response matrix
-    fullgamma <-if (nlevels(rhs[,nlevs]) == 1)
-        TRUE else FALSE
-#    if (!fullgamma && !global)
-#        warning("gamma diversity value might be meaningless")
-    ftmp <- vector("list", nlevs)
-    for (i in 1:nlevs) {
-        ftmp[[i]] <- as.formula(paste("~", tlab[i], "- 1"))
-    }
-
-    ## is there a method/burnin/thin in ... ?
-    method <- if (is.null(list(...)$method))
-        "r2dtable" else list(...)$method
-    burnin <- if (is.null(list(...)$burnin))
-        0 else list(...)$burnin
-    thin <- if (is.null(list(...)$thin))
-        1 else list(...)$thin
-
-    ## evaluate other arguments
-    index <- match.arg(index)
-    divfun <- switch(index,
-        "renyi" = function(x) renyi(x, scales=scales, hill = TRUE),
-        "tsallis" = function(x) tsallis(x, scales=scales, hill = TRUE))
-
-    ## cluster membership determination
-    nrhs <- rhs
-    nrhs <- sapply(nrhs, as.numeric)
-    idcl <- function(i) {
-        h <- nrhs[,i]
-        l <- nrhs[,(i-1)]
-        sapply(unique(l), function(i) h[l==i][1])
-    }
-    id <- lapply(2:nlevs, idcl)
-
-    ## this is the function passed to oecosimu
-    if (global) {
-        wdivfun <- function(x) {
-            if (fullgamma) {
-                tmp <- lapply(1:(nlevs-1), function(i) t(model.matrix(ftmp[[i]], rhs)) %*% x)
-                tmp[[nlevs]] <- matrix(colSums(x), nrow = 1, ncol = ncol(x))
-            } else {
-                tmp <- lapply(1:nlevs, function(i) t(model.matrix(ftmp[[i]], rhs)) %*% x)
-            }
-            raw <- sapply(1:nlevs, function(i) divfun(tmp[[i]]))
-            a <- sapply(raw, mean)
-            G <- a[nlevs]
-            b <- sapply(1:(nlevs-1), function(i) G / a[i])
-            if (relative)
-                b <- b / sapply(raw[1:(nlevs-1)], length)
-            c(a, b)
-        }
-    } else {
-        wdivfun <- function(x) {
-            if (fullgamma) {
-                tmp <- lapply(1:(nlevs-1), function(i) t(model.matrix(ftmp[[i]], rhs)) %*% x)
-                tmp[[nlevs]] <- matrix(colSums(x), nrow = 1, ncol = ncol(x))
-            } else {
-                tmp <- lapply(1:nlevs, function(i) t(model.matrix(ftmp[[i]], rhs)) %*% x)
-            }
-            a <- sapply(1:nlevs, function(i) divfun(tmp[[i]]))
-            am <- lapply(1:(nlevs-1), function(i) {
-                    sapply(1:length(unique(id[[i]])), function(ii) {
-                        mean(a[[i]][id[[i]]==ii])
-                    })
-                })
-            b <- lapply(1:(nlevs-1), function(i) a[[(i+1)]] / am[[i]])
-            bmax <- lapply(id, function(i) table(i))
-            if (relative)
-                b <- lapply(1:(nlevs-1), function(i) b[[i]] / bmax[[i]])
-            c(sapply(a, mean), sapply(b, mean))
-        }
-    }
-    if (nsimul > 0) {
-            sim <- oecosimu(lhs, wdivfun, method = method, nsimul=nsimul,
-                burnin=burnin, thin=thin)
-        } else {
-            sim <- wdivfun(lhs)
-            tmp <- rep(NA, length(sim))
-            sim <- list(statistic = sim,
-                oecosimu = list(z = tmp, pval = tmp, method = NA, statistic = sim))
-        }
-    nam <- c(paste("alpha", 1:(nlevs-1), sep="."), "gamma",
-        paste("beta", 1:(nlevs-1), sep="."))
-    names(sim$statistic) <- attr(sim$oecosimu$statistic, "names") <- nam
-    attr(sim, "call") <- match.call()
-    attr(sim, "index") <- index
-    attr(sim, "scales") <- scales
-    attr(sim, "global") <- TRUE
-    attr(sim, "n.levels") <- nlevs
-    attr(sim, "terms") <- tlab
-    attr(sim, "model") <- rhs
-    class(sim) <- c("multipart", "list")
+    ## run simulations
+    sim <- multipart.default(lhs, rhs, index = index, scales = scales,
+                             global = global, relative = relative,
+                             nsimul = nsimul, ...)
+    call <- match.call()
+    call[[1]] <- as.name("multipart")
+    attr(sim, "call") <- call
     sim
 }
diff --git a/R/oecosimu.R b/R/oecosimu.R
index abc3df3..1f14cbe 100644
--- a/R/oecosimu.R
+++ b/R/oecosimu.R
@@ -128,6 +128,7 @@
     ind$oecosimu <- list(z = z, means = means, pval = p, simulated=simind,
                          method=method,
                          statistic = indstat, alternative = alternative)
+    attr(ind, "call") <- match.call()
     class(ind) <- c("oecosimu", class(ind))
     ind
 }
diff --git a/R/permutest.betadisper.R b/R/permutest.betadisper.R
index 0cda88a..cd3a1ff 100644
--- a/R/permutest.betadisper.R
+++ b/R/permutest.betadisper.R
@@ -4,10 +4,10 @@
     t.statistic <- function(x, y) {
         m <- length(x)
         n <- length(y)
-        xbar <- mean(x) ## .Internal(mean(x))
-        ybar <- mean(y) ## .Internal(mean(y))
-        xvar <- var(x)  ## .Internal(cov(x, NULL, 1, FALSE))
-        yvar <- var(y)  ## .Internal(cov(y, NULL, 1, FALSE))
+        xbar <- mean(x)
+        ybar <- mean(y)
+        xvar <- var(x)
+        yvar <- var(y)
         pooled <- sqrt(((m-1)*xvar + (n-1)*yvar) / (m+n-2))
         (xbar - ybar) / (pooled * sqrt(1/m + 1/n))
     }
@@ -33,7 +33,6 @@
                         x$distances[x$group == z[2]])})
     }
     for(i in seq(along = res[-1])) {
-        ##perm <- permuted.index2(nobs, control = control)
         perm <- shuffle(nobs, control = control)
         perm.resid <- resids[perm]
         f <- qr.fitted(mod.Q, perm.resid)
diff --git a/R/plot.envfit.R b/R/plot.envfit.R
index 4622700..9a4b04c 100644
--- a/R/plot.envfit.R
+++ b/R/plot.envfit.R
@@ -1,13 +1,36 @@
 `plot.envfit` <-
-    function (x, choices = c(1, 2), arrow.mul, at = c(0, 0),
+    function (x, choices = c(1, 2), labels, arrow.mul, at = c(0, 0),
               axis = FALSE, p.max = NULL, col = "blue", bg, add = TRUE, ...)
 {
     formals(arrows) <- c(formals(arrows), alist(... = ))
+    ## get labels
+    labs <- list("v" = rownames(x$vectors$arrows),
+                 "f" = rownames(x$factors$centroids))
+    ## Change labels if user so wishes
+    if (!missing(labels)) {
+        ## input list of "vectors" and/or "factors"
+        if (is.list(labels)) {
+            if (!is.null(labs$v) && !is.null(labels$vectors))
+                labs$v <- labels$vectors
+            if (!is.null(labs$f) && !is.null(labels$factors))
+                labs$f <- labels$factors
+        } else {
+            ## input vector: either vectors or factors must be NULL,
+            ## and the existing set of labels is replaced
+            if (!is.null(labs$v) && !is.null(labs$f))
+                stop("needs a list with both 'vectors' and 'factors' labels")
+            if (!is.null(labs$v))
+                labs$v <- labels
+            else
+                labs$f <- labels
+        }
+    }
     vect <- NULL
     if (!is.null(p.max)) {
         if (!is.null(x$vectors)) {
             take <- x$vectors$pvals <= p.max
             x$vectors$arrows <- x$vectors$arrows[take, , drop = FALSE]
+            labs$v <- labs$v[take]
             x$vectors$r <- x$vectors$r[take]
             if (nrow(x$vectors$arrows) == 0)
                 x$vectors <- vect <- NULL
@@ -18,6 +41,7 @@
             take <- x$factors$var.id %in% nam
             x$factors$centroids <- x$factors$centroids[take,
                                                        , drop = FALSE]
+            labs$f <- labs$f[take]
             if (nrow(x$factors$centroids) == 0)
                 x$factors <- NULL
         }
@@ -42,16 +66,14 @@
         plot.new() ## needed for string widths and heights
         if(!is.null(vect)) {
             ## compute axis limits allowing space for labels
-            labs <- rownames(x$vectors$arrows)
-            sw <- strwidth(labs, ...)
-            sh <- strheight(labs, ...)
+            sw <- strwidth(labs$v, ...)
+            sh <- strheight(labs$v, ...)
             xlim <- range(at[1], vtext[,1] + sw, vtext[,1] - sw)
             ylim <- range(at[2], vtext[,2] + sh, vtext[,2] - sh)
             if(!is.null(x$factors)) {
                 ## if factors, also need to consider them
-                labs <- rownames(x$factors$centroids)
-                sw <- strwidth(labs, ...)
-                sh <- strheight(labs, ...)
+                sw <- strwidth(labs$f, ...)
+                sh <- strheight(labs$f, ...)
                 xlim <- range(xlim, x$factors$centroids[, choices[1]] + sw,
                               x$factors$centroids[, choices[1]] - sw)
                 ylim <- range(ylim, x$factors$centroids[, choices[2]] + sh,
@@ -66,9 +88,8 @@
             alabs <- colnames(vect)
             title(..., ylab = alabs[2], xlab = alabs[1])
         } else if (!is.null(x$factors)) {
-            labs <- rownames(x$factors$centroids)
-            sw <- strwidth(labs, ...)
-            sh <- strheight(labs, ...)
+            sw <- strwidth(labs$f, ...)
+            sh <- strheight(labs$f, ...)
             xlim <- range(at[1], x$factors$centroids[, choices[1]] + sw,
                           x$factors$centroids[, choices[1]] - sw)
             ylim <- range(at[2], x$factors$centroids[, choices[2]] + sh,
@@ -87,19 +108,17 @@
         arrows(at[1], at[2], vect[, 1], vect[, 2], len = 0.05,
                col = col)
         if (missing(bg))
-            text(vtext, rownames(x$vectors$arrows), col = col, ...)
+            text(vtext, labs$v, col = col, ...)
         else
-            ordilabel(vtext, labels = rownames(x$vectors$arrows),
-                      col = col, fill = bg, ...)
+            ordilabel(vtext, labels = labs$v, col = col, fill = bg, ...)
     }
     if (!is.null(x$factors)) {
         if (missing(bg))
             text(x$factors$centroids[, choices, drop = FALSE],
-                 rownames(x$factors$centroids), col = col, ...)
+                 labs$f, col = col, ...)
         else
             ordilabel(x$factors$centroids[, choices, drop = FALSE],
-                      labels = rownames(x$factors$centroids),
-                      col = col, fill = bg, ...)
+                      labels = labs$f, col = col, fill = bg, ...)
     }
     if (axis && !is.null(vect)) {
         axis(3, at = ax + at[1], labels = c(maxarr, 0, maxarr),
diff --git a/R/plot.radfit.R b/R/plot.radfit.R
index 501fb69..5149ec6 100644
--- a/R/plot.radfit.R
+++ b/R/plot.radfit.R
@@ -1,8 +1,13 @@
-"plot.radfit" <-
+`plot.radfit` <-
     function (x, BIC = FALSE, legend = TRUE, ...) 
 {
     if (length(x$y) == 0)
         stop("No species, nothing to plot")
+    ## if 'type = "n"', do not add legend (other types are not
+    ## supported)
+    type <- match.call(expand.dots = FALSE)$...$type
+    if (is.null(type))
+        type <- ""
     out <- plot(x$y, ...)
     if (length(x$y) == 1)
         return(invisible(out))
@@ -14,7 +19,7 @@
     lwd <- rep(1, ncol(fv))
     lwd[emph] <- 3
     matlines(fv, lty = 1, lwd = lwd, ...)
-    if (legend) {
+    if (legend && type != "n") {
         nm <- names(x$models)
         legend("topright", legend = nm, lty = 1, lwd = lwd, col = 1:6)
     }
diff --git a/R/points.radline.R b/R/points.radline.R
index 4861955..cf647d0 100644
--- a/R/points.radline.R
+++ b/R/points.radline.R
@@ -1,4 +1,4 @@
-"points.radline" <-
+`points.radline` <-
     function (x, ...) 
 {
     poi <- x$y
@@ -8,3 +8,9 @@
     class(out) <- "ordiplot"
     invisible(out)
 }
+
+`points.radfit` <-
+    function(x, ...)
+{
+    points.radline(x, ...)
+}
diff --git a/R/predict.radline.R b/R/predict.radline.R
new file mode 100644
index 0000000..80b1db1
--- /dev/null
+++ b/R/predict.radline.R
@@ -0,0 +1,50 @@
+### predict method for radline, radfit & radfit.frame
+
+### All functions take 'newdata' argument which need not be integer:
+### the functions can interpolate, but not necessarily extrapolate, or
+### the extrapolations may be NaN.
+
+`predict.radline`  <-
+    function(object, newdata, total, ...)
+{
+    ## newdata can be ranks
+    if (missing(newdata))
+        x <- seq_along(object$y)
+    else
+        x <- drop(as.matrix(newdata))
+    ## total number of individuals in the community
+    if (missing(total))
+        total <- sum(object$y)
+    ## adjustment for chagned total in call
+    adj <- total/sum(object$y)
+    nobs <- length(object$y)
+    p <- coef(object)
+    switch(object$model,
+           ## linear interpolation, no extrapolation
+           `Brokenstick` = approx(seq_len(nobs),
+           object$fitted.values, x, ...)$y * adj,
+           `Preemption` = exp(log(total) + log(p) + log(1 - p)*(x-1)),
+           ## NaN when rank outside proportional rank 0...1 
+           `Log-Normal` = {
+               slope <- diff(range(ppoints(nobs)))/(nobs-1)
+               intcpt <- 0.5 - slope * (nobs + 1) / 2
+               xnorm <- -qnorm(intcpt + slope * x)
+               exp(p[1] + p[2]*xnorm)*adj
+           },
+           `Zipf` = exp(log(total) + log(p[1]) + p[2]*log(x)),
+           `Zipf-Mandelbrot` = exp(log(total) + log(p[1]) +
+           p[2]*log(x + p[3]))
+           )
+}
+
+`predict.radfit`<-
+    function(object, ...)
+{
+    sapply(object$models, predict, ...)
+}
+
+`predict.radfit.frame`  <-
+    function(object, ...)
+{
+    lapply(object, predict, ...)
+}
diff --git a/R/print.adipart.R b/R/print.adipart.R
deleted file mode 100644
index 7278e56..0000000
--- a/R/print.adipart.R
+++ /dev/null
@@ -1,32 +0,0 @@
-print.adipart <-
-function(x, ...)
-{
-    n <- if (is.null(x$oecosimu$simulated))
-        0 else ncol(x$oecosimu$simulated)
-    if (n > 0)
-        cat("adipart with", n, "simulations using method",
-            dQuote(x$oecosimu$method), "\n")
-    else
-        cat("adipart ")
-    att <- attributes(x)
-    att$names <- att$call <- att$class <- att$n.levels <- att$terms <- att$model <- NULL
-    cat("with", paste(names(att), att, collapse=", "))
-
-    cat("\n\n")
-    cl <- class(x)
-    if (length(cl) > 1 && cl[2] != "list") {
-        NextMethod("print", x)
-        cat("\n")
-    }
-    if (!is.null(x$oecosimu$simulated)) {
-        tmp <- x$oecosimu$simulated
-    } else {
-        tmp <- data.matrix(x$oecosimu$statistic)
-    }
-    qu <- apply(tmp, 1, quantile, probs=c(0.025, 0.5, 0.975))
-    m <- cbind("statistic" = x$oecosimu$statistic,
-               "z" = x$oecosimu$z, t(qu),
-               "Pr(sim.)"=x$oecosimu$pval)
-    printCoefmat(m, ...)
-    invisible(x)
-}
diff --git a/R/print.hiersimu.R b/R/print.hiersimu.R
deleted file mode 100644
index 2c4cc17..0000000
--- a/R/print.hiersimu.R
+++ /dev/null
@@ -1,16 +0,0 @@
-print.hiersimu <-
-function (x, ...)
-{
-    cat("hiersimu with", ncol(x$oecosimu$simulated), "simulations\n\n")
-    cl <- class(x)
-    if (length(cl) > 1 && cl[2] != "list") {
-        NextMethod("print", x)
-        cat("\n")
-    }
-    qu <- apply(x$oecosimu$simulated, 1, quantile, probs = c(0.025,
-        0.5, 0.975))
-    m <- cbind(statistic = x$oecosimu$statistic, z = x$oecosimu$z,
-        t(qu), `Pr(sim.)` = x$oecosimu$pval)
-    printCoefmat(m, ...)
-    invisible(x)
-}
diff --git a/R/print.multipart.R b/R/print.multipart.R
deleted file mode 100644
index 68b49eb..0000000
--- a/R/print.multipart.R
+++ /dev/null
@@ -1,27 +0,0 @@
-print.multipart <-
-function(x, ...)
-{
-    n <- if (is.null(x$oecosimu$simulated))
-        0 else ncol(x$oecosimu$simulated)
-    cat("multipart with", n, "simulations\n")
-    att <- attributes(x)
-    att$names <- att$call <- att$class <- att$n.levels <- att$terms <- att$model <- NULL
-    cat("with", paste(names(att), att, collapse=", "))
-    cat("\n\n")
-    cl <- class(x)
-    if (length(cl) > 1 && cl[2] != "list") {
-        NextMethod("print", x)
-        cat("\n")
-    }
-    if (!is.null(x$oecosimu$simulated)) {
-        tmp <- x$oecosimu$simulated
-    } else {
-        tmp <- data.matrix(x$oecosimu$statistic)
-    }
-    qu <- apply(tmp, 1, quantile, probs=c(0.025, 0.5, 0.975))
-    m <- cbind("statistic" = x$oecosimu$statistic,
-               "z" = x$oecosimu$z, t(qu),
-               "Pr(sim.)"=x$oecosimu$pval)
-    printCoefmat(m, ...)
-    invisible(x)
-}
diff --git a/R/print.oecosimu.R b/R/print.oecosimu.R
index 48c1029..45b8100 100644
--- a/R/print.oecosimu.R
+++ b/R/print.oecosimu.R
@@ -1,12 +1,16 @@
 `print.oecosimu` <-
     function(x, ...)
 {
+    xx <- x ## return unmodified input object
     attr(x$oecosimu$method, "permfun") <- NULL
-    cat("oecosimu with", ncol(x$oecosimu$simulated), "simulations\n")
-    cat("simulation method", x$oecosimu$method)
+    cat(as.character(attr(x,"call")[[1]]), "object\n\n")
+    writeLines(strwrap(pasteCall(attr(x, "call"))))
+    cat("\n")
+    cat("simulation method", x$oecosimu$method, "with",
+        ncol(x$oecosimu$simulated), "simulations\n")
     if (length(att <- attributes(x$oecosimu$simulated)) > 1) {
         att$dim <- NULL
-        cat(" with", paste(names(att), att, collapse=", "))
+        cat("options: ", paste(names(att), att, collapse=", "))
     }
     alt.char <- switch(x$oecosimu$alternative,
                        two.sided = "not equal to",
@@ -17,9 +21,10 @@
 
     cat("\n\n")
     cl <- class(x)
-    if (length(cl) > 1 && cl[2] != "list") {
-        NextMethod("print", x)
-        cat("\n")
+    if ((length(cl) > 1 && cl[2] != "list" ) &&
+        !any(cl %in% c("adipart", "hiersimu", "multipart"))) {
+            NextMethod("print", x)
+            cat("\n")
     }
     probs <- switch(x$oecosimu$alternative,
                     two.sided = c(0.025, 0.5, 0.975),
@@ -35,5 +40,7 @@
         cat("\nNumber of NA cases removed from simulations:\n",
             nacount, "\n")
     }
-    invisible(x)   
+    invisible(xx)   
 }
+
+
diff --git a/R/protest.R b/R/protest.R
index ee6cba0..2ea5a53 100644
--- a/R/protest.R
+++ b/R/protest.R
@@ -12,9 +12,7 @@
         tmp <- procrustes(X, Y[take, ], symmetric = TRUE)$ss
         perm[i] <- sqrt(1 - tmp)
     }
-    perm <- c(sol$t0, perm)
-    permutations <- permutations + 1
-    Pval <- sum(perm >= sol$t0)/permutations
+    Pval <- (sum(perm >= sol$t0) + 1)/(permutations + 1)
     if (!missing(strata)) {
         strata <- deparse(substitute(strata))
         s.val <- strata
diff --git a/R/rad.null.R b/R/rad.null.R
index 9a187e1..d0b71dc 100644
--- a/R/rad.null.R
+++ b/R/rad.null.R
@@ -21,6 +21,7 @@
     }
     residuals <- x - fit
     rdf <- nsp
+    names(fit) <- names(x)
     p <- NA
     names(p) <- "S"
     out <- list(model = "Brokenstick", family=fam, y = x, coefficients = p,
diff --git a/R/rad.preempt.R b/R/rad.preempt.R
index 3b7f98d..3959ad4 100644
--- a/R/rad.preempt.R
+++ b/R/rad.preempt.R
@@ -45,6 +45,7 @@
             aic <- NA
         rdf <- length(x) - 1
     }
+    names(fit) <- names(x)
     names(p) <- c("alpha")
     out <- list(model = "Preemption", family = fam, y = x, coefficients = p, 
                 fitted.values = fit, aic = aic, rank = 1, df.residual = rdf, 
diff --git a/R/radfit.R b/R/radfit.R
index 98b88b1..c31e852 100644
--- a/R/radfit.R
+++ b/R/radfit.R
@@ -1,5 +1,5 @@
-"radfit" <-
-    function (...)
+`radfit` <-
+    function (x, ...)
 {
     UseMethod("radfit")
 }
diff --git a/R/radfit.data.frame.R b/R/radfit.data.frame.R
index c13b4e7..1c4427a 100644
--- a/R/radfit.data.frame.R
+++ b/R/radfit.data.frame.R
@@ -1,9 +1,9 @@
-"radfit.data.frame" <-
-    function(df, ...)
+`radfit.data.frame` <-
+    function(x, ...)
 {
-    ## df *must* have rownames
-    rownames(df) <- rownames(df, do.NULL = TRUE)
-    out <- apply(df, 1, radfit, ...)
+    ## x *must* have rownames
+    rownames(x) <- rownames(x, do.NULL = TRUE)
+    out <- apply(x, 1, radfit, ...)
     if (length(out) == 1)
         out <- out[[1]]
     else {
@@ -12,3 +12,9 @@
     }
     out
 }
+
+`radfit.matrix` <-
+    function(x, ...)
+{
+    radfit(as.data.frame(x), ...)
+}
diff --git a/R/radlattice.R b/R/radlattice.R
index acfd972..5c8c0bf 100644
--- a/R/radlattice.R
+++ b/R/radlattice.R
@@ -1,6 +1,8 @@
 `radlattice` <-
     function(x, BIC = FALSE, ...)
 {
+    if (!inherits(x, "radfit"))
+        stop("function only works with 'radfit' results for single site")
     require(lattice) || stop("requires package 'lattice'")
     y <- x$y
     fv <- unlist(fitted(x))
@@ -10,7 +12,11 @@
     Abundance <- rep(y, p)
     Rank <- rep(1:n, p)
     Model <- factor(rep(mods, each=n), levels = mods)
-    aic <- AIC(x, BIC = BIC)
+    if (BIC)
+        k <- log(length(y))
+    else
+        k <- 2
+    aic <- AIC(x, k = k)
     col <- trellis.par.get("superpose.line")$col
     if (length(col) > 1)
         col <- col[2]
diff --git a/inst/ChangeLog b/inst/ChangeLog
index 83fdb8c..88c5157 100644
--- a/inst/ChangeLog
+++ b/inst/ChangeLog
@@ -1,7 +1,41 @@
-$Date: 2012-06-17 21:05:31 +0300 (Sun, 17 Jun 2012) $
+$Date: 2012-10-07 19:50:53 +0300 (Sun, 07 Oct 2012) $
 
 VEGAN RELEASE VERSIONS at http://cran.r-project.org/
 
+Version 2.0-5 (released October 8, 2012)
+
+	* merge r2309: anova.cca.Rd edits.
+	* merge r2307: no line breaks within \code{} in Rd.
+	* merge r2305: proofread Rd files a..b.
+	* merge r2299: fix broken \link{}s in docs. 
+	* merge r2297: upgrade docs for L&L 2012 (3r ed.)
+	* merge r2291 thru 2296, 2298, 2300: radfit upgrade.
+	* merge r2287,8: scoping in anova.ccabyaxis and anova.ccabyterm.
+	* merge r2285: add predict.radfit.
+	* merge r2262, 2268:2270, and also r1950: mantel and
+	mantel.partial gained argument na.rm = FALSE. This needed hand
+	editing of merges, and also merging old r1950: beware and test.
+	* merge r2271: tweak varpart.Rd. (plot coloured circles).
+	* merge r2267: tweak vegdist.Rd (ref to vegdist).
+	* merge r2260: streamline adonis (internal changes).
+	* merge r2258: protect betadisper against changes in R-devel API.
+	* merge r2254: stylistic in examples of Rd files.
+	* merge r2252,56,57,61,64,66: add density methods for vegan
+	permutation functions and add plot.densityvegan to display these.
+	* merge r2250: do not use paste0 in envift.Rd (fails in R 2.14).
+	* merge r2249: fix vignette building with TeXLive 2012.
+	* merge r2246: remove dead code from cIndexKM() and R-devel CMD
+	check warning.
+	* merge r2244: more portable doc/Makefile
+	* merge r2237 thru 2240: add labels.envfit() and "labels" arg to
+	plot.envfit(). 
+	* merge r2227 thru 2235: adipart, hiersimu, multipart code
+	refactoring (especially formula method) and making to inherit from
+	oecosimu in printing the results.  The merge did not apply quite
+	cleanly, but oecosimu, print.oecosimu and NAMESPACE needed manual
+	editing (beware).
+	* merge r2225: biplot.rda 'type' fix.
+
 Version 2.0-4 (released June 18, 2012)
 
 	* merge r2215: plot.envfit() gains args 'bg' for background colour
diff --git a/inst/NEWS.Rd b/inst/NEWS.Rd
index adc48b5..427b7c3 100644
--- a/inst/NEWS.Rd
+++ b/inst/NEWS.Rd
@@ -2,6 +2,98 @@
 \title{vegan News}
 \encoding{UTF-8}
 
+\section{Changes in version 2.0-5}{
+
+  \subsection{BUG FIXES}{
+
+    \itemize{
+
+      \item \code{anova(<cca_object>, ...)} failed with \code{by =
+      "axis"} and \code{by = "term"}. The bug was reported by Dr Sven
+      Neulinger (Christian Albrecht University, Kiel, Germany).
+
+      \item \code{radlattice} did not honour argument \code{BIC = TRUE},
+      but always displayed AIC.
+
+    }
+
+  } % bug fixes
+
+  \subsection{NEW FUNCTIONS}{
+    \itemize{
+
+      \item Most \pkg{vegan} functions with permutation tests have now a
+      \code{density} method that can be used to find empirical
+      probability distributions of permutations. There is a new
+      \code{plot} method for these functions that displays both the
+      density and the observed statistic. The \code{density} function is
+      available for \code{adonis}, \code{anosim}, \code{mantel},
+      \code{mantel.partial}, \code{mrpp}, \code{permutest.cca} and
+      \code{procrustes}.
+
+      Function \code{adonis} can return several statistics, and it has
+      now a \code{densityplot} method (based on \pkg{lattice}).
+
+      Function \code{oecosimu} already had \code{density} and
+      \code{densityplot}, but they are now similar to other \pkg{vegan}
+      methods, and also work with \code{adipart}, \code{hiersimu} and
+      \code{multipart}.
+
+      \item \code{radfit} functions got a \code{predict} method that
+      also accepts arguments \code{newdata} and \code{total} for new
+      ranks and site totals for prediction.  The functions can also
+      interpolate to non-integer \dQuote{ranks}, and in some models
+      also extrapolate.
+
+    }
+  } % new functions
+
+  \subsection{NEW FEATURES}{
+    \itemize{
+
+      \item Labels can now be set in the \code{plot} of \code{envfit}
+      results. The labels must be given in the same order that the
+      function uses internally, and new support function \code{labels}
+      can be used to display the default labels in their correct order.
+
+      \item Mantel tests (functions \code{mantel} and
+      \code{mantel.partial}) gained argument \code{na.rm} which can be
+      used to remove missing values. This options should be used with
+      care: Permutation tests can be biased if the missing values were
+      originally in matching or fixed positions.
+
+      \item \code{radfit} results can be consistently accessed with
+      the same methods whether they were a single model for a single
+      site, all models for a single site or all models for all sites
+      in the data.  All functions now have methods \code{AIC},
+      \code{coef}, \code{deviance}, \code{logLik}, \code{fitted},
+      \code{predict} and \code{residuals}.
+
+    }
+
+  } % new features
+  
+  \subsection{INSTALLATION AND BUILDING}{
+    \itemize{
+
+      \item Building of \pkg{vegan} vignettes failed with the latest
+      version of LaTeX (TeXLive 2012).
+
+      \item \R{} versions later than 2.15-1 (including development
+      version) report warnings and errors when installing and checking
+      \pkg{vegan}, and you must upgrade \pkg{vegan} to this version.
+      The warnings concern functions \code{cIndexKM} and
+      \code{betadisper}, and the error occurs in \code{betadisper}.
+      These errors and warnings were triggered by internal changes in
+      \R.
+      
+    }
+  } % installation and building
+  
+
+} % version 2.0-5
+
+
 \section{Changes in version 2.0-4}{
 
   \subsection{BUG FIXES}{
diff --git a/inst/doc/FAQ-vegan.pdf b/inst/doc/FAQ-vegan.pdf
index 2534c79..985f70c 100644
Binary files a/inst/doc/FAQ-vegan.pdf and b/inst/doc/FAQ-vegan.pdf differ
diff --git a/inst/doc/Makefile b/inst/doc/Makefile
index e71e279..9e6c692 100644
--- a/inst/doc/Makefile
+++ b/inst/doc/Makefile
@@ -1,14 +1,14 @@
 all: FAQ-vegan.pdf decision-vegan.pdf intro-vegan.pdf diversity-vegan.pdf NEWS.html
 FAQ-vegan.pdf: FAQ-vegan.texi
-	texi2dvi --pdf --clean FAQ-vegan.texi
+	"$(R_HOME)/bin/R" CMD texi2dvi --pdf --clean FAQ-vegan.texi
 decision-vegan.pdf: decision-vegan.tex
-	texi2dvi --pdf --clean decision-vegan.tex
+	"$(R_HOME)/bin/R" CMD texi2dvi --pdf --clean decision-vegan.tex
 	-rm -f decision-vegan-0*.* Rplots.*
 intro-vegan.pdf: intro-vegan.tex
-	texi2dvi --pdf --clean intro-vegan.tex
+	"$(R_HOME)/bin/R" CMD texi2dvi --pdf --clean intro-vegan.tex
 	-rm -f intro-vegan-0*.* Rplots.*
 diversity-vegan.pdf: diversity-vegan.tex
-	texi2dvi --pdf --clean diversity-vegan.tex
+	"$(R_HOME)/bin/R" CMD texi2dvi --pdf --clean diversity-vegan.tex
 	-rm -f diversity-vegan-0*.* Rplots.*
 NEWS.html: ../NEWS.Rd
 	"$(R_HOME)/bin/R" CMD Rd2txt -t html ../NEWS.Rd -o NEWS.html
diff --git a/inst/doc/NEWS.html b/inst/doc/NEWS.html
index c4a17b7..c479a71 100644
--- a/inst/doc/NEWS.html
+++ b/inst/doc/NEWS.html
@@ -8,6 +8,115 @@
 
 <h2>vegan News</h2>
 
+<h3>Changes in version 2.0-5</h3>
+
+
+
+
+<h4>BUG FIXES</h4>
+
+
+
+<ul>
+<li> <p><code>anova(<cca_object>, ...)</code> failed with <code>by =
+      "axis"</code> and <code>by = "term"</code>. The bug was reported by Dr Sven
+Neulinger (Christian Albrecht University, Kiel, Germany).
+</p>
+</li>
+<li> <p><code>radlattice</code> did not honour argument <code>BIC = TRUE</code>,
+but always displayed AIC.
+</p>
+</li></ul>
+
+ 
+
+
+<h4>NEW FUNCTIONS</h4>
+
+
+
+<ul>
+<li><p> Most <span class="pkg">vegan</span> functions with permutation tests have now a
+<code>density</code> method that can be used to find empirical
+probability distributions of permutations. There is a new
+<code>plot</code> method for these functions that displays both the
+density and the observed statistic. The <code>density</code> function is
+available for <code>adonis</code>, <code>anosim</code>, <code>mantel</code>,
+<code>mantel.partial</code>, <code>mrpp</code>, <code>permutest.cca</code> and
+<code>procrustes</code>.
+</p>
+<p>Function <code>adonis</code> can return several statistics, and it has
+now a <code>densityplot</code> method (based on <span class="pkg">lattice</span>).
+</p>
+<p>Function <code>oecosimu</code> already had <code>density</code> and
+<code>densityplot</code>, but they are now similar to other <span class="pkg">vegan</span>
+methods, and also work with <code>adipart</code>, <code>hiersimu</code> and
+<code>multipart</code>.
+</p>
+</li>
+<li> <p><code>radfit</code> functions got a <code>predict</code> method that
+also accepts arguments <code>newdata</code> and <code>total</code> for new
+ranks and site totals for prediction.  The functions can also
+interpolate to non-integer “ranks”, and in some models
+also extrapolate.
+</p>
+</li></ul>
+
+ 
+
+
+<h4>NEW FEATURES</h4>
+
+
+
+<ul>
+<li><p> Labels can now be set in the <code>plot</code> of <code>envfit</code>
+results. The labels must be given in the same order that the
+function uses internally, and new support function <code>labels</code>
+can be used to display the default labels in their correct order.
+</p>
+</li>
+<li><p> Mantel tests (functions <code>mantel</code> and
+<code>mantel.partial</code>) gained argument <code>na.rm</code> which can be
+used to remove missing values. This options should be used with
+care: Permutation tests can be biased if the missing values were
+originally in matching or fixed positions.
+</p>
+</li>
+<li> <p><code>radfit</code> results can be consistently accessed with
+the same methods whether they were a single model for a single
+site, all models for a single site or all models for all sites
+in the data.  All functions now have methods <code>AIC</code>,
+<code>coef</code>, <code>deviance</code>, <code>logLik</code>, <code>fitted</code>,
+<code>predict</code> and <code>residuals</code>.
+</p>
+</li></ul>
+
+ 
+
+
+<h4>INSTALLATION AND BUILDING</h4>
+
+
+
+<ul>
+<li><p> Building of <span class="pkg">vegan</span> vignettes failed with the latest
+version of LaTeX (TeXLive 2012).
+</p>
+</li>
+<li> <p><font face="Courier New,Courier" color="#666666"><b>R</b></font> versions later than 2.15-1 (including development
+version) report warnings and errors when installing and checking
+<span class="pkg">vegan</span>, and you must upgrade <span class="pkg">vegan</span> to this version.
+The warnings concern functions <code>cIndexKM</code> and
+<code>betadisper</code>, and the error occurs in <code>betadisper</code>.
+These errors and warnings were triggered by internal changes in
+<font face="Courier New,Courier" color="#666666"><b>R</b></font>.
+</p>
+</li></ul>
+
+ 
+
+
 <h3>Changes in version 2.0-4</h3>
 
 
diff --git a/inst/doc/decision-vegan.pdf b/inst/doc/decision-vegan.pdf
index 099ab42..8844aa5 100644
Binary files a/inst/doc/decision-vegan.pdf and b/inst/doc/decision-vegan.pdf differ
diff --git a/inst/doc/decision-vegan.tex b/inst/doc/decision-vegan.tex
index 2d8d799..52ca7cc 100644
--- a/inst/doc/decision-vegan.tex
+++ b/inst/doc/decision-vegan.tex
@@ -23,8 +23,8 @@ another document.
 %% hijack Address for version info
 \Address{$ $Id: decision-vegan.Rnw 1709 2011-08-10 15:48:21Z jarioksa $ $
   processed with vegan
-2.0-4
-in R Under development (unstable) (2012-06-18 r59572) on \today}
+2.0-5
+in R Under development (unstable) (2012-10-08 r60901) on \today}
 \Footername{About this version}
 
 %% need no \usepackage{Sweave.sty}
@@ -500,17 +500,17 @@ Call: cca(formula = varespec[i, ] ~ Al + K, data = varechem)
 
               Inertia Proportion Rank
 Total         2.08320    1.00000     
-Constrained   0.19245    0.09238    2
-Unconstrained 1.89075    0.90762   21
+Constrained   0.13087    0.06282    2
+Unconstrained 1.95233    0.93718   21
 Inertia is mean squared contingency coefficient 
 
 Eigenvalues for constrained axes:
    CCA1    CCA2 
-0.10480 0.08764 
+0.07534 0.05553 
 
 Eigenvalues for unconstrained axes:
     CA1     CA2     CA3     CA4     CA5     CA6     CA7     CA8 
-0.46635 0.33696 0.23103 0.18183 0.14526 0.12137 0.09830 0.08624 
+0.49651 0.34903 0.23403 0.18364 0.17438 0.11899 0.10508 0.08591 
 (Showed only 8 of all 21 unconstrained eigenvalues)
 \end{Soutput}
 \end{Schunk}
@@ -537,11 +537,11 @@ R> proc <- procrustes(scores(tmp1, dis="lc", choi=1:14), scores(tmp2, dis="lc",
 R> max(residuals(proc))
 \end{Sinput}
 \begin{Soutput}
-[1] 2.703843e-14
+[1] 2.572443e-14
 \end{Soutput}
 \end{Schunk}
 In \code{cca} the difference would be somewhat larger than now
-observed 2.7038e-14 because site
+observed 2.5724e-14 because site
 weights used for environmental variables are shuffled with the species
 data.
 
diff --git a/inst/doc/diversity-vegan.Rnw b/inst/doc/diversity-vegan.Rnw
index b65c053..2195ec3 100644
--- a/inst/doc/diversity-vegan.Rnw
+++ b/inst/doc/diversity-vegan.Rnw
@@ -30,7 +30,7 @@
   pool}
 
 %% misuse next for scm data
-\Address{$ $Id: diversity-vegan.Rnw 1801 2011-09-07 16:04:10Z jarioksa $ $
+\Address{$ $Id: diversity-vegan.Rnw 2259 2012-08-23 15:13:58Z jarioksa $ $
   processed with vegan \Sexpr{packageDescription("vegan", field="Version")}
   in \Sexpr{R.version.string} on \today}
 \Footername{About this version}
@@ -49,9 +49,8 @@ options("prompt" = "R> ", "continue" = "+  ")
 
 \tableofcontents
 
-\section*{~}
-
-The \pkg{vegan} package has two major components:
+ ~\\[2ex]
+\noindent The \pkg{vegan} package has two major components:
 multivariate analysis (mainly ordination), and methods for diversity
 analysis of ecological communities.  This document gives an
 introduction to the latter.  Ordination methods are covered in other
diff --git a/inst/doc/diversity-vegan.pdf b/inst/doc/diversity-vegan.pdf
index 620a38e..8c4b1ed 100644
Binary files a/inst/doc/diversity-vegan.pdf and b/inst/doc/diversity-vegan.pdf differ
diff --git a/inst/doc/diversity-vegan.tex b/inst/doc/diversity-vegan.tex
index 4db4836..4678e6c 100644
--- a/inst/doc/diversity-vegan.tex
+++ b/inst/doc/diversity-vegan.tex
@@ -30,9 +30,9 @@
   pool}
 
 %% misuse next for scm data
-\Address{$ $Id: diversity-vegan.Rnw 1801 2011-09-07 16:04:10Z jarioksa $ $
-  processed with vegan 2.0-4
-  in R Under development (unstable) (2012-06-18 r59572) on \today}
+\Address{$ $Id: diversity-vegan.Rnw 2259 2012-08-23 15:13:58Z jarioksa $ $
+  processed with vegan 2.0-5
+  in R Under development (unstable) (2012-10-08 r60901) on \today}
 \Footername{About this version}
 
 %% need no \usepackage{Sweave}
@@ -42,9 +42,8 @@
 
 \tableofcontents
 
-\section*{~}
-
-The \pkg{vegan} package has two major components:
+ ~\\[2ex]
+\noindent The \pkg{vegan} package has two major components:
 multivariate analysis (mainly ordination), and methods for diversity
 analysis of ecological communities.  This document gives an
 introduction to the latter.  Ordination methods are covered in other
@@ -353,16 +352,16 @@ R> fish
 \end{Sinput}
 \begin{Soutput}
 Fisher log series model
-No. of species: 102 
+No. of species: 100 
 
       Estimate Std. Error
-alpha   44.064     5.5838
+alpha   40.996      5.183
 \end{Soutput}
 \end{Schunk}
 \begin{SCfigure}
 \includegraphics{diversity-vegan-018}
 \caption{Fisher's log-series fitted to one randomly selected site
-  (41).}
+  (20).}
 \label{fig:fisher}
 \end{SCfigure}
 We already saw $\alpha$ as a diversity index.  Now we also obtained
@@ -380,7 +379,7 @@ R> confint(fish)
 \end{Sinput}
 \begin{Soutput}
    2.5 %   97.5 % 
-34.21748 56.25208 
+31.84517 52.29146 
 \end{Soutput}
 \end{Schunk}
 
@@ -404,7 +403,7 @@ the lower octave.
 Function \code{prestondistr} directly
 maximizes truncated log-normal likelihood without binning data, and it
 is the recommended alternative.  Log-normal models  usually fit poorly
-to the BCI data, but here our random plot (number 41):
+to the BCI data, but here our random plot (number 20):
 \begin{Schunk}
 \begin{Sinput}
 R> prestondistr(BCI[k,])
@@ -412,15 +411,15 @@ R> prestondistr(BCI[k,])
 \begin{Soutput}
 Preston lognormal model
 Method: maximized likelihood to log2 abundances 
-No. of species: 102 
+No. of species: 100 
 
       mode      width         S0 
- 0.6954283  1.7260910 28.1654616 
+ 0.9849822  1.6642110 27.1323499 
 
 Frequencies by Octave
-                0        1        2       3        4        5
-Observed 22.50000 32.00000 21.00000 13.5000 8.000000 5.000000
-Fitted   25.96983 27.73039 21.16784 11.5513 4.506289 1.256727
+                0        1        2      3         4        5         6
+Observed 19.00000 29.00000 22.00000 12.500 13.500000 3.000000 1.0000000
+Fitted   22.77303 27.13125 22.52739 13.036  5.257389 1.477706 0.2894664
 \end{Soutput}
 \end{Schunk}
 
@@ -459,20 +458,20 @@ R> rad
 \end{Sinput}
 \begin{Soutput}
 RAD models, family poisson 
-No. of species 102, total abundance 402
-
-           par1      par2    par3    Deviance AIC      BIC     
-Null                                  59.9800 342.5388 342.5388
-Preemption  0.040337                  47.8801 332.4389 335.0639
-Lognormal   0.80833   1.0772          24.2781 310.8369 316.0869
-Zipf        0.11928  -0.7903          36.6726 323.2314 328.4813
-Mandelbrot  2.3922   -1.5273  9.0795   7.0992 295.6580 303.5330
+No. of species 100, total abundance 429
+
+           par1      par2     par3   Deviance AIC     BIC    
+Null                                  45.067  331.756 331.756
+Preemption  0.039794                  33.728  322.417 325.022
+Lognormal   0.93253   1.0384          19.018  309.707 314.917
+Zipf        0.11174  -0.76181         42.503  333.192 338.402
+Mandelbrot  8.5711   -1.7894   14.87   7.286  299.975 307.790
 \end{Soutput}
 \end{Schunk}
 \begin{SCfigure}
 \includegraphics{diversity-vegan-022}
 \caption{Ranked abundance distribution models for a random plot
-  (no. 41).  The best model has the lowest \textsc{aic}.}
+  (no. 20).  The best model has the lowest \textsc{aic}.}
 \label{fig:rad}
 \end{SCfigure}
 
@@ -777,10 +776,10 @@ R> s <- sample(nrow(BCI), 25)
 R> specpool(BCI[s,])
 \end{Sinput}
 \begin{Soutput}
-    Species     chao  chao.se  jack1 jack1.se    jack2     boot
-All     204 219.5588 8.479327 226.08 6.505136 232.2517 214.9225
-     boot.se  n
-All 3.946063 25
+    Species    chao  chao.se  jack1 jack1.se    jack2     boot  boot.se
+All     207 222.125 8.427135 228.12 7.533605 234.2533 217.6358 4.392338
+     n
+All 25
 \end{Soutput}
 \end{Schunk}
 
@@ -797,12 +796,12 @@ two of these methods:
 R> estimateR(BCI[k,])
 \end{Sinput}
 \begin{Soutput}
-                 41
-S.obs    102.000000
-S.chao1  151.500000
-se.chao1  21.334825
-S.ACE    161.646487
-se.ACE     6.839367
+                 20
+S.obs    100.000000
+S.chao1  133.476190
+se.chao1  15.441519
+S.ACE    147.592810
+se.ACE     6.621289
 \end{Soutput}
 \end{Schunk}
 Chao's method is similar as above, but uses another, ``unbiased''
@@ -845,14 +844,14 @@ R> veiledspec(prestondistr(BCI[k,]))
 \end{Sinput}
 \begin{Soutput}
 Extrapolated     Observed       Veiled 
-   121.86261    102.00000     19.86261 
+   113.18418    100.00000     13.18418 
 \end{Soutput}
 \begin{Sinput}
 R> veiledspec(BCI[k,])
 \end{Sinput}
 \begin{Soutput}
 Extrapolated     Observed       Veiled 
-   146.08199    102.00000     44.08199 
+   122.24639    100.00000     22.24639 
 \end{Soutput}
 \end{Schunk}
 
diff --git a/inst/doc/intro-vegan.Rnw b/inst/doc/intro-vegan.Rnw
index 15f3641..143689a 100644
--- a/inst/doc/intro-vegan.Rnw
+++ b/inst/doc/intro-vegan.Rnw
@@ -34,7 +34,7 @@
   vector, fitted environmental surface, permutation tests}
 
 %% misuse of the address field for revision data
-\Address{$ $Id: intro-vegan.Rnw 1709 2011-08-10 15:48:21Z jarioksa $ $
+\Address{$ $Id: intro-vegan.Rnw 2259 2012-08-23 15:13:58Z jarioksa $ $
   processed with vegan
 \Sexpr{packageDescription("vegan", field="Version")}
 in \Sexpr{R.version.string} on \today}
@@ -55,9 +55,8 @@ options("prompt" = "R> ", "continue" = "+  ")
 
 \tableofcontents
 
-\section*{~}
-
-\pkg{Vegan} is a package for community ecologists.  This
+~\\[2ex]
+\noindent \pkg{Vegan} is a package for community ecologists.  This
 documents explains how the commonly used ordination methods can be
 performed in \pkg{vegan}.  The document only is a very basic
 introduction.  Another document (\emph{vegan tutorial})
diff --git a/inst/doc/intro-vegan.pdf b/inst/doc/intro-vegan.pdf
index a30de5e..634a4aa 100644
Binary files a/inst/doc/intro-vegan.pdf and b/inst/doc/intro-vegan.pdf differ
diff --git a/inst/doc/intro-vegan.tex b/inst/doc/intro-vegan.tex
index f1b8862..dd03cae 100644
--- a/inst/doc/intro-vegan.tex
+++ b/inst/doc/intro-vegan.tex
@@ -34,10 +34,10 @@
   vector, fitted environmental surface, permutation tests}
 
 %% misuse of the address field for revision data
-\Address{$ $Id: intro-vegan.Rnw 1709 2011-08-10 15:48:21Z jarioksa $ $
+\Address{$ $Id: intro-vegan.Rnw 2259 2012-08-23 15:13:58Z jarioksa $ $
   processed with vegan
-2.0-4
-in R Under development (unstable) (2012-06-18 r59572) on \today}
+2.0-5
+in R Under development (unstable) (2012-10-08 r60901) on \today}
 \Footername{About this version}
 
 %% need no \usepackage{Sweave}
@@ -48,9 +48,8 @@ in R Under development (unstable) (2012-06-18 r59572) on \today}
 
 \tableofcontents
 
-\section*{~}
-
-\pkg{Vegan} is a package for community ecologists.  This
+~\\[2ex]
+\noindent \pkg{Vegan} is a package for community ecologists.  This
 documents explains how the commonly used ordination methods can be
 performed in \pkg{vegan}.  The document only is a very basic
 introduction.  Another document (\emph{vegan tutorial})
@@ -128,8 +127,15 @@ R> ord <- metaMDS(dune)
 \end{Sinput}
 \begin{Soutput}
 Run 0 stress 0.1192691 
-Run 1 stress 0.11927 
-... procrustes: rmse 0.001094789  max resid 0.003155725 
+Run 1 stress 0.1812939 
+Run 2 stress 0.1808915 
+Run 3 stress 0.1809588 
+Run 4 stress 0.1183188 
+... New best solution
+... procrustes: rmse 0.02021248  max resid 0.06425558 
+Run 5 stress 0.1192686 
+Run 6 stress 0.1183202 
+... procrustes: rmse 0.0005415232  max resid 0.001655342 
 *** Solution reached
 \end{Soutput}
 \begin{Sinput}
@@ -145,9 +151,9 @@ Data:     dune
 Distance: bray 
 
 Dimensions: 2 
-Stress:     0.1192691 
+Stress:     0.1183188 
 Stress type 1, weak ties
-Two convergent solutions found after 1 tries
+Two convergent solutions found after 6 tries
 Scaling: centring, PC rotation, halfchange scaling 
 Species: expanded scores based on ‘dune’ 
 \end{Soutput}
@@ -305,8 +311,8 @@ R> ord.fit
 \begin{Soutput}
 ***VECTORS
 
-     NMDS1   NMDS2   r2  Pr(>r)  
-A1 0.99055 0.13712 0.38 0.01598 *
+     NMDS1   NMDS2    r2  Pr(>r)  
+A1 0.96467 0.26345 0.365 0.02398 *
 ---
 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
 P values based on 1000 permutations.
@@ -315,14 +321,14 @@ P values based on 1000 permutations.
 
 Centroids:
                NMDS1   NMDS2
-ManagementBF -0.4478 -0.0190
-ManagementHF -0.2695 -0.1255
-ManagementNM  0.2995  0.5797
-ManagementSF  0.1490 -0.4656
+ManagementBF -0.4534 -0.0103
+ManagementHF -0.2638 -0.1282
+ManagementNM  0.2957  0.5790
+ManagementSF  0.1509 -0.4670
 
 Goodness of fit:
                r2   Pr(>r)   
-Management 0.4142 0.008991 **
+Management 0.4134 0.008991 **
 ---
 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
 P values based on 1000 permutations.
@@ -350,12 +356,12 @@ Link function: identity
 
 Formula:
 y ~ s(x1, x2, k = knots)
-<environment: 0x4501508>
+<environment: 0x3525cd8>
 
 Estimated degrees of freedom:
 2  total = 3 
 
-GCV score: 3.871584
+GCV score: 3.964875
 \end{Soutput}
 \end{Schunk}
 \begin{SCfigure}
@@ -496,7 +502,7 @@ Terms added sequentially (first to last)
 
 Model: cca(formula = dune ~ A1 + Management, data = dune.env)
            Df  Chisq      F N.Perm Pr(>F)   
-A1          1 0.2248 2.5245    199  0.020 * 
+A1          1 0.2248 2.5245    199  0.015 * 
 Management  3 0.5550 2.0780    199  0.005 **
 Residual   15 1.3355                        
 ---
@@ -517,7 +523,7 @@ Marginal effects of terms
 
 Model: cca(formula = dune ~ A1 + Management, data = dune.env)
            Df  Chisq      F N.Perm  Pr(>F)   
-A1          1 0.1759 1.9761   1199 0.03333 * 
+A1          1 0.1759 1.9761    699 0.02571 * 
 Management  3 0.5550 2.0780    199 0.00500 **
 Residual   15 1.3355                         
 ---
@@ -534,9 +540,9 @@ R> anova(ord, by="axis", perm=500)
 Model: cca(formula = dune ~ A1 + Management, data = dune.env)
          Df  Chisq      F N.Perm Pr(>F)   
 CCA1      1 0.3187 3.5801    199  0.005 **
-CCA2      1 0.2372 2.6640    199  0.010 **
-CCA3      1 0.1322 1.4845    299  0.100 . 
-CCA4      1 0.0917 1.0297     99  0.380   
+CCA2      1 0.2372 2.6640    199  0.015 * 
+CCA3      1 0.1322 1.4845     99  0.160   
+CCA4      1 0.0917 1.0297     99  0.390   
 Residual 15 1.3355                        
 ---
 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
@@ -591,10 +597,10 @@ Permutation test for cca under reduced model
 Terms added sequentially (first to last)
 
 Model: cca(formula = dune ~ A1 + Management + Condition(Moisture), data = dune.env)
-           Df  Chisq      F N.Perm Pr(>F)   
-A1          1 0.1154 1.4190     99   0.10 . 
-Management  3 0.3954 1.6205     99   0.01 **
-Residual   12 0.9761                        
+           Df  Chisq      F N.Perm Pr(>F)  
+A1          1 0.1154 1.4190     99   0.10 .
+Management  3 0.3954 1.6205     99   0.03 *
+Residual   12 0.9761                       
 ---
 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
 \end{Soutput}
@@ -613,7 +619,7 @@ Permutations stratified within 'Moisture'
 
 Model: cca(formula = dune ~ A1 + Management + Condition(Moisture), data = dune.env)
            Df  Chisq      F N.Perm Pr(>F)   
-A1          1 0.1154 1.4190     99   0.26   
+A1          1 0.1154 1.4190     99   0.27   
 Management  3 0.3954 1.6205     99   0.01 **
 Residual   12 0.9761                        
 ---
diff --git a/man/BCI.Rd b/man/BCI.Rd
index 35bce5d..fd51327 100644
--- a/man/BCI.Rd
+++ b/man/BCI.Rd
@@ -21,7 +21,7 @@
   The data frame contains only the Barro Colorado Island subset of the
   original data.
 
-  The quadrats are located in a regular grid. See\code{examples} for the
+  The quadrats are located in a regular grid. See \code{examples} for the
   coordinates. 
 
 }
diff --git a/man/RsquareAdj.Rd b/man/RsquareAdj.Rd
index 86ef42b..0cd90f7 100644
--- a/man/RsquareAdj.Rd
+++ b/man/RsquareAdj.Rd
@@ -48,8 +48,9 @@ Adjusted R-square
   \eqn{R^2}{R-squared}. The adjusted \eqn{R^2}{R-squared} is found as
   the difference of adjusted \eqn{R^2}{R-squared} values of joint effect
   of partial and constraining terms and partial term alone, and it is
-  the same as the adjusted \eqn{R^2}{R-squared} of component \code{[a] =
-  X1|X2} in two-component variation partition in \code{\link{varpart}}.
+  the same as the adjusted \eqn{R^2}{R-squared} of component 
+  \code{[a] = X1|X2} in two-component variation partition in 
+  \code{\link{varpart}}.
   }
 
 \value{ The functions return a list of items \code{r.squared} and
diff --git a/man/add1.cca.Rd b/man/add1.cca.Rd
index 4ba8d13..c4167d6 100644
--- a/man/add1.cca.Rd
+++ b/man/add1.cca.Rd
@@ -4,7 +4,7 @@
 
 \title{Add or Drop  Single Terms to a Constrained Ordination Model }
 \description{
-Compute all single terms that can be added or dropped from a
+Compute all single terms that can be added to or dropped from a
 constrained ordination model.
 }
 \usage{
@@ -19,7 +19,7 @@ constrained ordination model.
   \code{\link{cca}}, \code{\link{rda}} or \code{\link{capscale}}. }
   \item{scope}{ A formula giving the terms to be considered for adding
   or dropping; see \code{\link{add1}} for details.}
-  \item{test}{ Should a permutation test added using \code{\link{anova.cca}}. }
+  \item{test}{ Should a permutation test be added using \code{\link{anova.cca}}. }
   \item{pstep}{Number of permutations in one step, passed as argument
   \code{step} to \code{\link{anova.cca}}.}
   \item{perm.max}{ Maximum number of permutation in \code{\link{anova.cca}}. }
@@ -36,16 +36,16 @@ constrained ordination model.
   Function \code{add1.cca} will implement a test for single term
   additions that is not directly available in \code{\link{anova.cca}}.
 
-  Functions are used implicitly in \code{\link{step}} and
-  \code{\link{ordistep}}. The \code{\link{deviance.cca}} and
-  \code{\link{deviance.rda}} used in \code{\link{step}} have no firm
-  basis, and setting argument \code{test = "permutation"} may help in
-  getting useful insight into validity of model building. Function
-  \code{\link{ordistep}} calls alternately \code{drop1.cca} and
-  \code{add1.cca} with argument \code{test = "permutation"} and
-  selects variables by their permutation \eqn{P}-values.  Meticulous
-  use of \code{add1.cca} and \code{drop1.cca} will allow more
-  judicious model building.
+  Functions are used implicitly in \code{\link{step}},
+  \code{\link{ordiR2step}} and \code{\link{ordistep}}. The
+  \code{\link{deviance.cca}} and \code{\link{deviance.rda}} used in
+  \code{\link{step}} have no firm basis, and setting argument \code{test
+  = "permutation"} may help in getting useful insight into validity of
+  model building. Function \code{\link{ordistep}} calls alternately
+  \code{drop1.cca} and \code{add1.cca} with argument 
+  \code{test = "permutation"} and selects variables by their permutation
+  \eqn{P}-values.  Meticulous use of \code{add1.cca} and
+  \code{drop1.cca} will allow more judicious model building.
 
   The default \code{perm.max} is set to a low value, because
   permutation tests can take a long time. It should be sufficient to
@@ -62,7 +62,7 @@ constrained ordination model.
 
 \seealso{ \code{\link{add1}}, \code{\link{drop1}} and
   \code{\link{anova.cca}} for basic methods. You probably need these
-  functions with \code{\link{step}} and \code{link{ordistep}}. Functions
+  functions with \code{\link{step}} and \code{\link{ordistep}}. Functions
   \code{\link{deviance.cca}} and \code{\link{extractAIC.cca}} are used
   to produce the other arguments than test results in the
   output. Functions \code{\link{cca}}, \code{\link{rda}} and
diff --git a/man/adipart.Rd b/man/adipart.Rd
index 3a0f775..851b321 100644
--- a/man/adipart.Rd
+++ b/man/adipart.Rd
@@ -3,11 +3,10 @@
 \alias{adipart}
 \alias{adipart.default}
 \alias{adipart.formula}
-\alias{print.adipart}
 \alias{hiersimu}
 \alias{hiersimu.default}
 \alias{hiersimu.formula}
-\alias{print.hiersimu}
+
 \title{Additive Diversity Partitioning and Hierarchical Null Model Testing}
 \description{
 In additive diversity partitioning, mean values of alpha diversity at lower levels of a sampling 
@@ -33,89 +32,122 @@ hiersimu(...)
   \item{x}{A matrix with same number of rows as in \code{y}, columns
     coding the levels of sampling hierarchy. The number of groups within
     the hierarchy must decrease from left to right. If \code{x} is missing,
-    two levels are assumed: each row is a group in the first level, and
-    all rows are in the same group in the second level.}
-  \item{formula}{A two sided model formula in the form \code{y ~ x}, where \code{y} 
-    is the community data matrix with samples as rows and species as column. Right 
-    hand side (\code{x}) must contain factors referring to levels of sampling hierarchy, 
-    terms from right to left will be treated as nested (first column is the lowest, 
-    last is the highest level). These variables must be factors in order to unambiguous 
-    handling. Interaction terms are not allowed.}
-  \item{data}{A data frame where to look for variables defined in the right hand side 
-    of \code{formula}. If missing, variables are looked in the global environment.}
+    function performs an overall decomposition into alpha, beta and
+    gamma diversities.}
+  \item{formula}{A two sided model formula in the form \code{y ~ x},
+    where \code{y} is the community data matrix with samples as rows and
+    species as column. Right hand side (\code{x}) must grouping vaiables
+    referring to levels of sampling hierarchy, terms from right to left
+    will be treated as nested (first column is the lowest, last is the
+    highest level, at least two levels specified). Interaction terms are
+    not allowed.}
+
+  \item{data}{A data frame where to look for variables defined in the
+    right hand side of \code{formula}. If missing, variables are looked
+    in the global environment.}
+
   \item{index}{Character, the diversity index to be calculated (see Details).}
-  \item{weights}{Character, \code{"unif"} for uniform weights, \code{"prop"} for 
-    weighting proportional to sample abundances to use in weighted averaging of individual 
-    alpha values within strata of a given level of the sampling hierarchy.}
-  \item{relative}{Logical, if \code{TRUE} then alpha and beta diversity values are given 
-    relative to the value of gamma for function \code{adipart}.}
-  \item{nsimul}{Number of permutation to use if \code{matr} is not of class 'permat'.
-    If \code{nsimul = 0}, only the \code{FUN} argument is evaluated. It is thus possible
-    to reuse the statistic values without using a null model.}
-  \item{FUN}{A function to be used by \code{hiersimu}. This must be fully specified,
-    because currently other arguments cannot be passed to this function via \code{\dots}.}
-  \item{location}{Character, identifies which function (mean or median) is to be used to 
-    calculate location of the samples.}
-  \item{drop.highest}{Logical, to drop the highest level or not. When \code{FUN} 
-    evaluates only arrays with at least 2 dimensions, highest level should be dropped, 
-    or not selected at all.}
-  \item{\dots}{Other arguments passed to functions, e.g. base of logarithm for 
-    Shannon diversity, or \code{method}, \code{thin} or \code{burnin} arguments for
-    \code{\link{oecosimu}}.}
+
+  \item{weights}{Character, \code{"unif"} for uniform weights,
+    \code{"prop"} for weighting proportional to sample abundances to use
+    in weighted averaging of individual alpha values within strata of a
+    given level of the sampling hierarchy.}
+
+  \item{relative}{Logical, if \code{TRUE} then alpha and beta diversity
+    values are given relative to the value of gamma for function
+    \code{adipart}.}
+
+  \item{nsimul}{Number of permutations to use if \code{matr} is not of
+    class 'permat'.  If \code{nsimul = 0}, only the \code{FUN} argument
+    is evaluated. It is thus possible to reuse the statistic values
+    without using a null model.}
+
+  \item{FUN}{A function to be used by \code{hiersimu}. This must be
+    fully specified, because currently other arguments cannot be passed
+    to this function via \code{\dots}.}
+
+  \item{location}{Character, identifies which function (mean or median)
+    is to be used to calculate location of the samples.}
+
+  \item{drop.highest}{Logical, to drop the highest level or not. When
+    \code{FUN} evaluates only arrays with at least 2 dimensions, highest
+    level should be dropped, or not selected at all.}
+
+  \item{\dots}{Other arguments passed to functions, e.g. base of
+    logarithm for Shannon diversity, or \code{method}, \code{thin} or
+    \code{burnin} arguments for \code{\link{oecosimu}}.}
 }
+
 \details{
-Additive diversity partitioning means that mean alpha and beta diversity adds up to gamma 
-diversity, thus beta diversity is measured in the same dimensions as alpha and gamma 
-(Lande 1996). This additive procedure is than extended across multiple scales in a 
-hierarchical sampling design with \eqn{i = 1, 2, 3, \ldots, m} levels of sampling 
-(Crist et al. 2003). Samples in lower hierarchical levels are nested within higher level 
-units, thus from \eqn{i=1} to \eqn{i=m} grain size is increasing under constant survey 
-extent. At each level \eqn{i}, \eqn{\alpha_i} denotes average diversity found within samples.
-
-At the highest sampling level, the diversity components are calculated as 
-\deqn{\beta_m = \gamma - \alpha_m}{beta_m = gamma - alpha_m} 
-For each lower sampling level as
-\deqn{\beta_i = \alpha_{i+1} - \alpha_i}{beta_i = alpha_i+1 - alpha_i}
-Then, the additive partition of diversity is 
-\deqn{\gamma = \alpha_1 + \sum_{i=1}^m \beta_i}{gamma = alpha_1 + sum(beta_i)}
-
-Average alpha components can be weighted uniformly (\code{weight="unif"}) to calculate 
-it as simple average, or proportionally to sample abundances (\code{weight="prop"}) to 
-calculate it as weighted average as follows
-\deqn{\alpha_i = \sum_{j=1}^{n_i} D_{ij} w_{ij}}{alpha_i = sum(D_ij*w_ij)}
-where \eqn{D_{ij}} is the diversity index and \eqn{w_{ij}} is the weight calculated for 
-the \eqn{j}th sample at the \eqn{i}th sampling level.
-
-The implementation of additive diversity partitioning in \code{adipart} follows Crist et 
-al. 2003. It is based on species richness (\eqn{S}, not \eqn{S-1}), Shannon's and 
-Simpson's diversity indices stated as the \code{index} argument.
-
-The expected diversity components are calculated \code{nsimul} times by individual based 
-randomisation of the community data matrix. This is done by the \code{"r2dtable"} method
-in \code{\link{oecosimu}} by default.
-
-\code{hiersimu} works almost the same as \code{adipart}, but without comparing the actual 
-statistic values returned by \code{FUN} to the highest possible value (cf. gamma diversity). 
-This is so, because in most of the cases, it is difficult to ensure additive properties of 
-the mean statistic values along the hierarchy.
+
+  Additive diversity partitioning means that mean alpha and beta
+  diversities add up to gamma diversity, thus beta diversity is measured
+  in the same dimensions as alpha and gamma (Lande 1996). This additive
+  procedure is then extended across multiple scales in a hierarchical
+  sampling design with \eqn{i = 1, 2, 3, \ldots, m} levels of sampling
+  (Crist et al. 2003). Samples in lower hierarchical levels are nested
+  within higher level units, thus from \eqn{i=1} to \eqn{i=m} grain size
+  is increasing under constant survey extent. At each level \eqn{i},
+  \eqn{\alpha_i} denotes average diversity found within samples.
+
+  At the highest sampling level, the diversity components are calculated
+  as \deqn{\beta_m  = \gamma -  \alpha_m}{beta_m = gamma -  alpha_m} For
+  each  lower   sampling  level   as  \deqn{\beta_i  =   \alpha_{i+1}  -
+  \alpha_i}{beta_i =  alpha_i+1 - alpha_i} Then,  the additive partition
+  of diversity is \deqn{\gamma  = \alpha_1 + \sum_{i=1}^m \beta_i}{gamma
+  = alpha_1 + sum(beta_i)}
+
+  Average alpha components can be weighted uniformly
+  (\code{weight="unif"}) to calculate it as simple average, or
+  proportionally to sample abundances (\code{weight="prop"}) to
+  calculate it as weighted average as follows \deqn{\alpha_i =
+  \sum_{j=1}^{n_i} D_{ij} w_{ij}}{alpha_i = sum(D_ij*w_ij)} where
+  \eqn{D_{ij}} is the diversity index and \eqn{w_{ij}} is the weight
+  calculated for the \eqn{j}th sample at the \eqn{i}th sampling level.
+
+  The implementation of additive diversity partitioning in
+  \code{adipart} follows Crist et al. 2003. It is based on species
+  richness (\eqn{S}, not \eqn{S-1}), Shannon's and Simpson's diversity
+  indices stated as the \code{index} argument.
+
+  The expected diversity components are calculated \code{nsimul} times
+  by individual based randomisation of the community data matrix. This
+  is done by the \code{"r2dtable"} method in \code{\link{oecosimu}} by
+  default.
+
+  \code{hiersimu} works almost in the same way as \code{adipart}, but
+  without comparing the actual statistic values returned by \code{FUN}
+  to the highest possible value (cf. gamma diversity).  This is so,
+  because in most of the cases, it is difficult to ensure additive
+  properties of the mean statistic values along the hierarchy.
+
 }
 \value{
-An object of class 'adipart' or 'hiersimu' with same structure as 'oecosimu' objects.
+
+  An object of class \code{"adipart"} or \code{"hiersimu"} with same
+  structure as \code{\link{oecosimu}} objects.
+
 }
+
 \references{
-Crist, T.O., Veech, J.A., Gering, J.C. and Summerville,
-K.S. (2003). Partitioning species diversity across landscapes and regions:
-a hierarchical analysis of \eqn{\alpha}, \eqn{\beta}, and
-\eqn{\gamma}-diversity.
-\emph{Am. Nat.}, \bold{162}, 734--743.
-
-Lande, R. (1996). Statistics and partitioning of species
-diversity, and similarity among multiple communities.
-\emph{Oikos}, \bold{76}, 5--13.
+
+  Crist,   T.O.,   Veech,    J.A.,   Gering,   J.C.   and   Summerville,
+  K.S.  (2003).  Partitioning species  diversity  across landscapes  and
+  regions:  a hierarchical  analysis of  \eqn{\alpha},  \eqn{\beta}, and
+  \eqn{\gamma}-diversity.  \emph{Am. Nat.}, \bold{162}, 734--743.
+
+  Lande, R.  (1996). Statistics and partitioning of species diversity,
+  and similarity among multiple communities.  \emph{Oikos}, \bold{76},
+  5--13.
+
 }
 
-\author{\enc{Péter Sólymos}{Peter Solymos}, \email{solymos at ualberta.ca}}
-\seealso{See \code{\link{oecosimu}} for permutation settings and calculating \eqn{p}-values.}
+\author{
+
+  \enc{Péter Sólymos}{Peter Solymos}, \email{solymos at ualberta.ca}}
+  \seealso{See \code{\link{oecosimu}} for permutation settings and
+  calculating \eqn{p}-values.  }
+
 \examples{
 ## NOTE: 'nsimul' argument usually needs to be >= 99
 ## here much lower value is used for demonstration
@@ -128,10 +160,10 @@ cutter <- function (x, cut = seq(0, 10, by = 2.5)) {
     out <- rep(1, length(x))
     for (i in 2:(length(cut) - 1))
         out[which(x > cut[i] & x <= cut[(i + 1)])] <- i
-    return(as.factor(out))}
+    return(out)}
 ## The hierarchy of sample aggregation
 levsm <- data.frame(
-    l1=as.factor(1:nrow(mite)),
+    l1=1:nrow(mite),
     l2=cutter(mite.xy$y, cut = seq(0, 10, by = 2.5)),
     l3=cutter(mite.xy$y, cut = seq(0, 10, by = 5)),
     l4=cutter(mite.xy$y, cut = seq(0, 10, by = 10)))
diff --git a/man/adonis.Rd b/man/adonis.Rd
index 899497f..b7a21c1 100644
--- a/man/adonis.Rd
+++ b/man/adonis.Rd
@@ -7,7 +7,7 @@
 \description{Analysis of variance using distance matrices --- for
   partitioning distance matrices among sources of variation and fitting
   linear models (e.g., factors, polynomial regression) to distance 
-  matrices; uses a permutation test with pseudo-F ratios.}
+  matrices; uses a permutation test with pseudo-\eqn{F} ratios.}
 
 \usage{
 adonis(formula, data, permutations = 999, method = "bray",
@@ -60,7 +60,7 @@ variance, Excoffier, Smouse, and Quattro, 1992;
 and nested factors.
 
 If the experimental design has nestedness, then use \code{strata} to
-test hypotheses. For instance, imagine we are testing the whether a
+test hypotheses. For instance, imagine we are testing whether a
 plant community is influenced by nitrate amendments, and we have two
 replicate plots at each of two levels of nitrate (0, 10 ppm). We have
 replicated the experiment in three fields with (perhaps) different
@@ -90,8 +90,8 @@ residuals. Permutations of the raw data may have better small sample
 characteristics. Further, the precise meaning of hypothesis tests will
 depend upon precisely what is permuted. The strata argument keeps groups
 intact for a particular hypothesis test where one does not want to
-permute the data among particular groups. For instance, \code{strata =
-B} causes permutations among levels of \code{A} but retains data within
+permute the data among particular groups. For instance, \code{strata = B} 
+causes permutations among levels of \code{A} but retains data within
 levels of \code{B} (no permutation among levels of \code{B}). See
 \code{\link{permutations}} for additional details on permutation tests
 in Vegan.
@@ -116,8 +116,8 @@ overview of rules.
 
   \item{aov.tab}{Typical AOV table showing sources of variation,
     degrees of freedom, sequential sums of squares, mean squares,
-    \eqn{F} statistics, partial R-squared and \eqn{P} values, based on \eqn{N}
-    permutations. }
+    \eqn{F} statistics, partial \eqn{R^2}{R-squared} and \eqn{P}
+    values, based on \eqn{N} permutations. }
   \item{coefficients}{ matrix of coefficients of the linear model, with
     rows representing sources of variation and columns representing
     species; each column represents a fit of a species abundance to the
@@ -129,10 +129,12 @@ overview of rules.
     sites; each column represents a fit of a sites distances (from all
     other sites) to the  linear model.These are what you get when you
     fit distances of one site to
-    your predictors. } 
+    your predictors. }   
   \item{f.perms}{ an \eqn{N} by \eqn{m} matrix of the null \eqn{F}
-    statistics for each source of variation  based on \eqn{N}
-    permutations of the data.}
+    statistics for each source of variation based on \eqn{N}
+    permutations of the data. The distribution of a single term can be
+    inspected with \code{\link{density.adonis}} function, or all terms
+    simultaneously with \code{densityplot.adonis}.}
   \item{model.matrix}{The \code{\link{model.matrix}} for the right hand
     side of the formula.}
   \item{terms}{The \code{\link{terms}} component of the model.}
@@ -143,7 +145,7 @@ overview of rules.
   by different within-group variation (dispersion) instead of different
   mean values of the groups (see Warton et al. 2012 for a general
   analysis). However, it seems that \code{adonis} is less sensitive to
-  dispersion effects than some of its alternatives (\code{link{anosim}},
+  dispersion effects than some of its alternatives (\code{\link{anosim}},
   \code{\link{mrpp}}). Function \code{\link{betadisper}} is a sister
   function to \code{adonis} to study the differences in dispersion
   within the same geometric framework.
@@ -211,10 +213,10 @@ ordihull(mod, group=dat$NO3, show="10", col=3)
 ordispider(mod, group=dat$field, lty=3, col="red")
 
 ### Correct hypothesis test (with strata)
-adonis(Y ~ NO3, data=dat, strata=dat$field, perm=1e3)
+adonis(Y ~ NO3, data=dat, strata=dat$field, perm=999)
 
 ### Incorrect (no strata)
-adonis(Y ~ NO3, data=dat, perm=1e3)
+adonis(Y ~ NO3, data=dat, perm=999)
 }
 
 \keyword{multivariate }
diff --git a/man/anosim.Rd b/man/anosim.Rd
index 3fb0217..390df72 100644
--- a/man/anosim.Rd
+++ b/man/anosim.Rd
@@ -21,7 +21,7 @@ anosim(dat, grouping, permutations = 999, distance = "bray", strata)
   \item{permutations}{Number of permutation to assess the significance
     of the ANOSIM statistic. }
   \item{distance}{Choice of distance metric that measures the
-    dissimilarity between two observations . See \code{\link{vegdist}} for
+    dissimilarity between two observations. See \code{\link{vegdist}} for
     options.  This will be used if \code{dat} was not a dissimilarity
     structure or a symmetric square matrix.}  
   \item{strata}{An integer vector or factor specifying the strata for
@@ -51,10 +51,11 @@ anosim(dat, grouping, permutations = 999, distance = "bray", strata)
   grouping.
 
   The statistical significance of observed \eqn{R} is assessed by
-  permuting the grouping vector to obtain the empirical
-  distribution of \eqn{R} under null-model.  See
-  \code{\link{permutations}} for additional details on permutation tests
-  in Vegan.
+  permuting the grouping vector to obtain the empirical distribution
+  of \eqn{R} under null-model.  See \code{\link{permutations}} for
+  additional details on permutation tests in Vegan. The distribution
+  of simulated values can be inspected with the \code{density}
+  function.
 
   The function has \code{summary} and \code{plot} methods.  These both
   show valuable information to assess the validity of the method:  The
@@ -69,7 +70,8 @@ anosim(dat, grouping, permutations = 999, distance = "bray", strata)
   \item{call }{Function call.}
   \item{statistic}{The value of ANOSIM statistic \eqn{R}}
   \item{signif}{Significance from permutation.}
-  \item{perm}{Permutation values of \eqn{R}}
+  \item{perm}{Permutation values of \eqn{R}. The distribution of
+    permutation values can be inspected with function \code{\link{density.anosim}}.}
   \item{class.vec}{Factor with value \code{Between} for dissimilarities
     between classes and class name for corresponding dissimilarity
     within class.}
diff --git a/man/anova.cca.Rd b/man/anova.cca.Rd
index 353d6b9..626bfaa 100644
--- a/man/anova.cca.Rd
+++ b/man/anova.cca.Rd
@@ -14,8 +14,8 @@
 \description{
   The function performs an ANOVA like permutation test for Constrained
   Correspondence Analysis (\code{\link{cca}}), Redundancy Analysis
-  (\code{\link{rda}}) or Constrained Analysis of Principal Coordinates
-  (\code{\link{capscale}}) to assess the significance of constraints.
+  (\code{\link{rda}}) or distance-based Redundancy Analysis
+  (dbRDA, \code{\link{capscale}}) to assess the significance of constraints.
 }
 \usage{
 \method{anova}{cca}(object, alpha=0.05, beta=0.01, step=100, perm.max=9999,
@@ -63,42 +63,6 @@ permutest(x, ...)
   Function \code{permutest.cca} is the proper workhorse, but
   \code{anova.cca} passes all parameters to \code{permutest.cca}.
 
-  In \code{anova.cca} the number of permutations is controlled by
-  targeted \dQuote{critical} \eqn{P} value (\code{alpha}) and accepted Type
-  II or rejection error (\code{beta}).  If the results of permutations
-  differ from the targeted \code{alpha} at risk level given by
-  \code{beta}, the permutations are
-  terminated.  If the current estimate of \eqn{P} does not
-  differ significantly from \code{alpha} of the alternative hypothesis,
-  the permutations are
-  continued with \code{step} new permutations (at the first step, the
-  number of permutations is \code{step - 1}).  However, with \code{by =
-    "terms"} a fixed number of permutations will be used, and this 
-  is given by argument \code{permutations}, or if this is missing,
-  by \code{step}.  
-  
-  The function \code{permutest.cca} implements a permutation test for
-  the \dQuote{significance} of constraints in \code{\link{cca}},
-  \code{\link{rda}} or \code{\link{capscale}}.  Community data are
-  permuted with choice \code{model = "direct"}, residuals after
-  partial CCA/RDA/CAP with choice \code{model = "reduced"} (default),
-  and residuals after CCA/RDA/CAP under choice \code{model = "full"}.
-  If there is no partial CCA/RDA/CAP stage, \code{model = "reduced"}
-  simply permutes the data and is equivalent to \code{model = "direct"}. 
-  The test statistic is ``pseudo-\eqn{F}'',
-  which is the ratio of constrained and unconstrained total Inertia
-  (Chi-squares, variances or something similar), each divided by their
-  respective ranks.  If there are no conditions (\dQuote{partial} terms), the
-  sum of all eigenvalues remains constant, so that pseudo-\eqn{F} and
-  eigenvalues would give equal results.  In partial CCA/RDA/CAP, the
-  effect of conditioning variables (\dQuote{covariables}) is removed before
-  permutation, and these residuals are added to the non-permuted fitted
-  values of partial CCA (fitted values of \code{X ~ Z}).  Consequently,
-  the total Chi-square is not fixed, and test based on pseudo-\eqn{F}
-  would differ from the test based on plain eigenvalues. CCA is a
-  weighted method, and environmental data are re-weighted at each
-  permutation step using permuted weights. 
-
   The default test is for the sum of all constrained eigenvalues.
   Setting \code{first = TRUE} will perform a test for the first
   constrained eigenvalue.  Argument \code{first} can be set either in
@@ -129,29 +93,65 @@ permutest(x, ...)
   will start from the same \code{\link{.Random.seed}}, and the seed
   will be advanced to the value after the longest permutation at the
   exit from the function.  
+  
+  In \code{anova.cca} the number of permutations is controlled by
+  targeted \dQuote{critical} \eqn{P} value (\code{alpha}) and accepted
+  Type II or rejection error (\code{beta}).  If the results of
+  permutations differ from the targeted \code{alpha} at risk level given
+  by \code{beta}, the permutations are terminated.  If the current
+  estimate of \eqn{P} does not differ significantly from \code{alpha} of
+  the alternative hypothesis, the permutations are continued with
+  \code{step} new permutations (at the first step, the number of
+  permutations is \code{step - 1}).  However, with \code{by="terms"} a
+  fixed number of permutations will be used, and this is given by
+  argument \code{permutations}, or if this is missing, by \code{step}.
+  
+  Community data are permuted with choice \code{model="direct"},
+  residuals after partial CCA/ RDA/ dbRDA with choice \code{model="reduced"} 
+  (default), and residuals after CCA/ RDA/ dbRDA under choice
+  \code{model="full"}.  If there is no partial CCA/ RDA/ dbRDA stage,
+  \code{model="reduced"} simply permutes the data and is equivalent to
+  \code{model="direct"}.  The test statistic is \dQuote{pseudo-\eqn{F}},
+  which is the ratio of constrained and unconstrained total Inertia
+  (Chi-squares, variances or something similar), each divided by their
+  respective ranks.  If there are no conditions (\dQuote{partial}
+  terms), the sum of all eigenvalues remains constant, so that
+  pseudo-\eqn{F} and eigenvalues would give equal results.  In partial
+  CCA/ RDA/ dbRDA, the effect of conditioning variables
+  (\dQuote{covariables}) is removed before permutation, and these
+  residuals are added to the non-permuted fitted values of partial CCA
+  (fitted values of \code{X ~ Z}).  Consequently, the total Chi-square
+  is not fixed, and test based on pseudo-\eqn{F} would differ from the
+  test based on plain eigenvalues. CCA is a weighted method, and
+  environmental data are re-weighted at each permutation step using
+  permuted weights.
 }
 
-\value{
+\value{ 
   Function \code{permutest.cca} returns an object of class
-  \code{"permutest.cca"}, which has its own \code{print} method.  The
-  function \code{anova.cca} calls \code{permutest.cca}, fills an
-  \code{\link{anova}} table and uses \code{\link{print.anova}} for printing.
+  \code{"permutest.cca"}, which has its own \code{print} method. The
+  distribution of permuted \eqn{F} values can be inspected with
+  \code{\link{density.permutest.cca}} function.  The function
+  \code{anova.cca} calls \code{permutest.cca}, fills an
+  \code{\link{anova}} table and uses \code{\link{print.anova}} for
+  printing.
 }
+
 \note{
   Some cases of \code{anova} need access to the original data on
   constraints (at least \code{by = "term"} and \code{by = "margin"}),
   and they may fail if data are unavailable.
  
   The default permutation \code{model} changed from \code{"direct"} to
-  \code{"reduced"} in \pkg{vegan} version 1.14-11 (release version
-  1.15-0), and you must explicitly set \code{model = "direct"} for
-  compatibility with the old version.
+  \code{"reduced"} in \pkg{vegan} version 1.15-0, and you must
+  explicitly set \code{model = "direct"} for compatibility with the old
+  version.
 
   Tests \code{by = "terms"} and \code{by = "margin"} are consistent
   only when \code{model = "direct"}.  
 }
 \references{
-  Legendre, P. and Legendre, L. (1998). \emph{Numerical Ecology}. 2nd
+  Legendre, P. and Legendre, L. (2012). \emph{Numerical Ecology}. 3rd
   English ed. Elsevier.
 
   Legendre, P., Oksanen, J. and ter Braak, C.J.F. (2011). Testing the
diff --git a/man/beals.Rd b/man/beals.Rd
index 4925e69..3fa0969 100644
--- a/man/beals.Rd
+++ b/man/beals.Rd
@@ -63,18 +63,17 @@ swan(x, maxit = Inf)
   conditioned probabilities and weighted averages.
 
   Beals smoothing was originally suggested as a method of data
-  transformation to remove excessive zeros (Beals 1984, McCune
-  1994).  However, it is not a suitable method for this purpose since it
-  does not maintain the information on species presences: a species may
-  have a higher probability of occurrence at a site where it does not
-  occur than at sites where it occurs. Moreover, it regularizes data
-  too strongly. The method may be useful in identifying species that
-  belong to the species pool (Ewald 2002) or to identify suitable
-  unoccupied patches in metapopulation analysis
-  (\enc{Münzbergová}{Munzbergova} & Herben 
-  2004). In this case, the function should be called with \code{include
-  = FALSE} for cross-validation smoothing for species; argument
-  \code{species} can be used if only one species is studied.
+  transformation to remove excessive zeros (Beals 1984, McCune 1994).
+  However, it is not a suitable method for this purpose since it does
+  not maintain the information on species presences: a species may have
+  a higher probability of occurrence at a site where it does not occur
+  than at sites where it occurs. Moreover, it regularizes data too
+  strongly. The method may be useful in identifying species that belong
+  to the species pool (Ewald 2002) or to identify suitable unoccupied
+  patches in metapopulation analysis (\enc{Münzbergová}{Munzbergova} &
+  Herben 2004). In this case, the function should be called with
+  \code{include=FALSE} for cross-validation smoothing for species;
+  argument \code{species} can be used if only one species is studied.
 
   Swan (1970) suggested replacing zero values with degrees of absence of
   a species in a community data matrix. Swan expressed the method in
diff --git a/man/betadisper.Rd b/man/betadisper.Rd
index b6fd06a..b34d087 100644
--- a/man/betadisper.Rd
+++ b/man/betadisper.Rd
@@ -49,7 +49,7 @@ betadisper(d, group, type = c("median","centroid"), bias.adjust = FALSE)
   \item{group}{vector describing the group structure, usually a factor
     or an object that can be coerced to a factor using
     \code{\link[base]{as.factor}}. Can consist of a factor with a single
-    level (i.e.~one group).}
+    level (i.e., one group).}
   \item{type}{the type of analysis to perform. Use the spatial median or
     the group centroid? The spatial median is now the default.}
   \item{bias.adjust}{logical: adjust for small sample bias in beta diversity estimates?}
@@ -100,10 +100,11 @@ betadisper(d, group, type = c("median","centroid"), bias.adjust = FALSE)
   \sqrt{\Delta^2(u_{ij}^+, c_i^+) - \Delta^2(u_{ij}^-, c_i^-)},}{z[ij]^c
   = sqrt(Delta^2(u[ij]^+, c[i]^+) - Delta^2(u[ij]^-, c[i]^-)),} where
   \eqn{\Delta^2}{Delta^2} is the squared Euclidean distance between
-  \eqn{u_{ij}}{u[ij]}, the principal coordinate for the \eqn{j^{th}}{jth}
-  point in the \eqn{i^{th}}{ith} group, and \eqn{c_i}{c[i]}, the
-  coordinate of the centroid for the \eqn{i^{th}}{ith} group. The
-  super-scripted \eqn{+} and \eqn{-} indicate the real and imaginary
+  \eqn{u_{ij}}{u[ij]}, the principal coordinate for the \eqn{j}th
+  point in the \eqn{i}th group, and \eqn{c_i}{c[i]}, the
+  coordinate of the centroid for the \eqn{i}th group. The
+  super-scripted \sQuote{\eqn{+}} and \sQuote{\eqn{-}} indicate
+  the real and imaginary
   parts respectively. This is equation (3) in Anderson (2006). If the
   imaginary part is greater in magnitude than the real part, then we
   would be taking the square root of a negative value, resulting in
@@ -114,9 +115,9 @@ betadisper(d, group, type = c("median","centroid"), bias.adjust = FALSE)
   
   To test if one or more groups is more variable than the others, ANOVA
   of the distances to group centroids can be performed and parametric
-  theory used to interpret the significance of F. An alternative is to
+  theory used to interpret the significance of \eqn{F}. An alternative is to
   use a permutation test. \code{\link{permutest.betadisper}} permutes model
-  residuals to generate a permutation distribution of F under the Null
+  residuals to generate a permutation distribution of \eqn{F} under the Null
   hypothesis of no difference in dispersion between groups.
 
   Pairwise comparisons of group mean dispersions can also be performed
@@ -186,12 +187,12 @@ betadisper(d, group, type = c("median","centroid"), bias.adjust = FALSE)
   prior to performing the analysis.
 }
 \section{Warning}{
-  Stewart Schultz noticed that the permutation test for \code{type =
-  "centroid"} had the wrong type I error and was anti-conservative. As
-  such, the default for \code{type} has been changed to \code{"median"},
-  which uses the spatial median as the group centroid. Tests suggests
-  that the permutation test for this type of analysis gives the correct
-  error rates.
+  Stewart Schultz noticed that the permutation test for
+  \code{type="centroid"} had the wrong type I error and was
+  anti-conservative. As such, the default for \code{type} has been
+  changed to \code{"median"}, which uses the spatial median as the group
+  centroid. Tests suggests that the permutation test for this type of
+  analysis gives the correct error rates.
 }
 \references{
   Anderson, M. J. (2001) A new method for non-parametric multivariate 
diff --git a/man/betadiver.Rd b/man/betadiver.Rd
index 892a8d5..c5fe882 100644
--- a/man/betadiver.Rd
+++ b/man/betadiver.Rd
@@ -52,11 +52,10 @@ betadiver(x, method = NA, order = FALSE, help = FALSE, ...)
 
   Function \code{betadiver} finds all indices reviewed by Koleff et
   al. (2003). All these indices could be found with function
-  \code{\link{designdist}} which uses different notation, but the
-  current function provides a conventional shortcut. The function
-  only finds the indices. The proper analysis must be done with
-  functions such as \code{\link{betadisper}}, \code{\link{adonis}} or
-  \code{\link{mantel}}. 
+  \code{\link{designdist}}, but the current function provides a
+  conventional shortcut. The function only finds the indices. The proper
+  analysis must be done with functions such as \code{\link{betadisper}},
+  \code{\link{adonis}} or \code{\link{mantel}}.
 
   The indices are directly taken from Table 1 of Koleff et al. (2003),
   and they can be selected either by the index number or the subscript
@@ -92,11 +91,6 @@ betadiver(x, method = NA, order = FALSE, help = FALSE, ...)
   are two such similarity indices.
 }
 
-\note{The argument \code{method} was called \code{index} in older
-  versions of the function (upto \pkg{vegan} version
-  1.17-11). Argument \code{index} is deprecated, but still recognized
-  with a warning. }
-
 \references{
 
   Baselga, A. (2010) Partitioning the turnover and nestedness
diff --git a/man/bgdispersal.Rd b/man/bgdispersal.Rd
index 0f80c0d..d1de766 100644
--- a/man/bgdispersal.Rd
+++ b/man/bgdispersal.Rd
@@ -6,7 +6,7 @@
 
 \description{ This function computes coefficients of dispersal direction
 between geographically connected areas, as defined by Legendre and
-Legendre (1984), and also described in Legendre and Legendre (1998,
+Legendre (1984), and also described in Legendre and Legendre (2012,
 section 13.3.4). }
 
 \usage{
@@ -48,22 +48,27 @@ and DD2 only are sought. }
 
 \value{
 Function \code{bgdispersal} returns a list containing the following matrices:
-\item{ DD1 }{ \eqn{DD1[j,k] = (a * (b - c))/((a + b + c)^2)} }
-\item{ DD2 }{ \eqn{DD2[j,k] = (2*a * (b - c))/((2*a + b + c) * (a + b +
+\item{ DD1 }{ \eqn{DD1_{j,k} = (a(b - c))/((a + b + c)^2)}{DD1[j,k] = (a * (b - c))/((a + b + c)^2)} }
+\item{ DD2 }{ \eqn{DD2_{j,k} = (2 a (b - c))/((2a + b + c)  (a + b +
+    c))}{DD2[j,k] = (2*a * (b - c))/((2*a + b + c) * (a + b +
     c))}
   where \eqn{a}, \eqn{b}, and \eqn{c} have the 
 same meaning as in the computation of binary 
 similarity coefficients. }
-\item{ DD3 }{ DD3[j,k] = \eqn{W*(A-B) / ((A+B-W)^2)} }
-\item{ DD4 }{ DD4[j,k] = \eqn{2*W*(A-B) / ((A+B)*(A+B-W)})
+\item{ DD3 }{ \eqn{DD3_{j,k} = {W(A-B) / (A+B-W)^2} }{DD3[j,k] = W*(A-B) / (A+B-W)^2} }
+\item{ DD4 }{ \eqn{DD4_{j,k} = 2W(A-B) / ((A+B)(A+B-W))}{DD4[j,k] = 2*W*(A-B) / ((A+B)*(A+B-W))}
 where \code{W = sum(pmin(vector1, vector2))}, \code{A = sum(vector1)},
 \code{B = sum(vector2)} }
+
 \item{ McNemar }{ McNemar chi-square statistic of asymmetry (Sokal and
-Rohlf 1995): \eqn{2*(b*log(b) + c*log(c) - (b+c)*log((b+c)/2)) / q}
-where \eqn{q = 1 + 1/(2*(b+c))} (Williams correction for continuity) }
+  Rohlf 1995):
+  \eqn{2(b \log(b) + c \log(c) - (b+c) \log((b+c)/2)) / q}{2*(b*log(b) + c*log(c) - (b+c)*log((b+c)/2)) / q},
+  where \eqn{q = 1 + 1/(2(b+c))}{q = 1 + 1/(2*(b+c))}
+  (Williams correction for continuity) }
 \item{ prob.McNemar }{ probabilities associated 
 with McNemar statistics, chi-square test. H0: no 
 asymmetry in \eqn{(b-c)}. }
+
 }
 
 
@@ -72,7 +77,7 @@ asymmetry in \eqn{(b-c)}. }
   freshwater fishes in the Québec
   peninsula. \emph{Can. J. Fish. Aquat. Sci.} \strong{41}: 1781-1802.
   
-  Legendre, P. and L. Legendre. 1998. \emph{Numerical ecology}, 2nd
+  Legendre, P. and L. Legendre. 2012. \emph{Numerical ecology}, 3rd
   English edition. Elsevier Science BV, Amsterdam.
   
   Sokal, R. R. and F. J. Rohlf. 1995. \emph{Biometry. The principles and
diff --git a/man/bioenv.Rd b/man/bioenv.Rd
index a23b1d7..faefe00 100644
--- a/man/bioenv.Rd
+++ b/man/bioenv.Rd
@@ -26,7 +26,7 @@
     in \code{\link{vegdist}}. This is ignored if \code{comm} are dissimilarities.}
   \item{upto}{Maximum number of parameters in studied subsets.}
   \item{formula, data}{Model \code{\link{formula}} and data.}
-  \item{trace}{Trace the advance of calculations }
+  \item{trace}{Trace the calculations }
   \item{partial}{Dissimilarities partialled out when inspecting
     variables in \code{env}.}
   \item{...}{Other arguments passed to \code{\link{cor}}.}
diff --git a/man/capscale.Rd b/man/capscale.Rd
index f7ecdcc..d18562d 100644
--- a/man/capscale.Rd
+++ b/man/capscale.Rd
@@ -56,7 +56,7 @@ capscale(formula, data, distance = "euclidean", sqrt.dist = FALSE,
      that all eigenvalues are non-negative in the underlying
      Principal Co-ordinates Analysis (see \code{\link{cmdscale}} 
      for details). This implements \dQuote{correction method 2} of
-     Legendre & Legendre (1998, p. 434). The negative eigenvalues are
+     Legendre & Legendre (2012, p. 503). The negative eigenvalues are
      caused by using semi-metric or non-metric dissimilarities with
      basically metric \code{\link{cmdscale}}. They are harmless and
      ignored in \code{capscale}, but you also can avoid warnings with
@@ -164,7 +164,7 @@ capscale(formula, data, distance = "euclidean", sqrt.dist = FALSE,
   analysis: testing multispecies responses in multifactorial ecological
   experiments. \emph{Ecological Monographs} 69, 1--24.
 
-  Legendre, P. & Legendre, L. (1998).  \emph{Numerical Ecology}. 2nd English
+  Legendre, P. & Legendre, L. (2012).  \emph{Numerical Ecology}. 3rd English
   Edition. Elsevier
 }
 \author{ Jari Oksanen }
diff --git a/man/cascadeKM.Rd b/man/cascadeKM.Rd
index e75e9b8..68a5404 100644
--- a/man/cascadeKM.Rd
+++ b/man/cascadeKM.Rd
@@ -117,7 +117,7 @@ cIndexKM(y, x, index = "all")
   distance matrix. (3) The first principal coordinate is used as the new
   order of the objects in the graph. A simplified algorithm is used to
   compute the first principal coordinate only, using the iterative
-  algorithm described in Legendre & Legendre (1998, Table 9.10). The
+  algorithm described in Legendre & Legendre (2012). The
   full distance matrix among objects is never computed; this avoids
   the problem of storing it when the number of objects is
   large. Distance values are computed as they are needed by the
@@ -157,7 +157,7 @@ cIndexKM(y, x, index = "all")
   methods used in multivariate analysis. \emph{Biometrika} \strong{53}:
   325-338.
   
-  Legendre, P. & L. Legendre. 1998. \emph{Numerical ecology}, 2nd
+  Legendre, P. & L. Legendre. 2012. \emph{Numerical ecology}, 3rd
   English edition. Elsevier Science BV, Amsterdam.
   
   Milligan, G. W. & M. C. Cooper. 1985. An examination of procedures for
diff --git a/man/cca.Rd b/man/cca.Rd
index 78bf59f..5d36232 100644
--- a/man/cca.Rd
+++ b/man/cca.Rd
@@ -60,7 +60,7 @@
   Functions \code{cca} and \code{rda} are  similar to popular
   proprietary software \code{Canoco}, although the implementation is
   completely different.  The functions are based on Legendre &
-  Legendre's (1998) algorithm: in \code{cca}
+  Legendre's (2012) algorithm: in \code{cca}
   Chi-square transformed data matrix is subjected to weighted linear
   regression on constraining variables, and the fitted values are
   submitted to correspondence analysis performed via singular value
@@ -170,7 +170,7 @@
 \references{ The original method was by ter Braak, but the current
   implementations follows Legendre and Legendre.
 
-  Legendre, P. and Legendre, L. (1998) \emph{Numerical Ecology}. 2nd English
+  Legendre, P. and Legendre, L. (2012) \emph{Numerical Ecology}. 3rd English
   ed. Elsevier.
 
   McCune, B. (1997) Influence of noisy environmental data on canonical
diff --git a/man/cca.object.Rd b/man/cca.object.Rd
index 212477e..9292e18 100644
--- a/man/cca.object.Rd
+++ b/man/cca.object.Rd
@@ -192,7 +192,7 @@
   \code{cca.object} as well). 
 }
 \references{
-  Legendre, P. and Legendre, L. (1998) \emph{Numerical Ecology}. 2nd English
+  Legendre, P. and Legendre, L. (2012) \emph{Numerical Ecology}. 3rd English
   ed. Elsevier.
 }
 \author{ Jari Oksanen }
diff --git a/man/clamtest.Rd b/man/clamtest.Rd
index d4a7d2c..0fcb282 100644
--- a/man/clamtest.Rd
+++ b/man/clamtest.Rd
@@ -141,9 +141,9 @@ The original method (Chazdon et al. 2011) has two major problems:
 \examples{
 data(mite)
 data(mite.env)
-x <- clamtest(mite, mite.env$Shrub=="None", alpha=0.005)
-summary(x)
-head(x)
-plot(x)
+sol <- clamtest(mite, mite.env$Shrub=="None", alpha=0.005)
+summary(sol)
+head(sol)
+plot(sol)
 }
 \keyword{ htest }
diff --git a/man/decostand.Rd b/man/decostand.Rd
index f543802..1d520ee 100644
--- a/man/decostand.Rd
+++ b/man/decostand.Rd
@@ -50,8 +50,7 @@ wisconsin(x)
     distance, the distances should be similar to the
     Chi-square distance used in correspondence analysis. However, the
     results from \code{\link{cmdscale}} would still differ, since
-    CA is a weighted ordination method (default \code{MARGIN =
-      1}).
+    CA is a weighted ordination method (default \code{MARGIN = 1}).
     \item \code{hellinger}: square root of \code{method = "total"}
     (Legendre & Gallagher 2001).
     \item \code{log}: logarithmic transformation as suggested by
@@ -61,9 +60,10 @@ wisconsin(x)
      to presences, and \code{logbase = Inf} gives the presence/absence
      scaling. Please note this is \emph{not} \eqn{\log(x+1)}{log(x+1)}.
      Anderson et al. (2006) suggested this for their (strongly) modified
-     Gower distance, but the standardization can be used independently
-     of distance indices.
-   }
+     Gower distance (implemented as \code{method = "altGower"} in 
+     \code{\link{vegdist}}), but the standardization can be used 
+     independently of distance indices. 
+  }
   Standardization, as contrasted to transformation, means that the
   entries are transformed relative to other entries.
 
diff --git a/man/density.adonis.Rd b/man/density.adonis.Rd
new file mode 100644
index 0000000..e6ec5f0
--- /dev/null
+++ b/man/density.adonis.Rd
@@ -0,0 +1,115 @@
+\name{density.adonis}
+\alias{density.adonis}
+\alias{density.anosim}
+\alias{density.mantel}
+\alias{density.mrpp}
+\alias{density.permutest.cca}
+\alias{density.protest}
+\alias{plot.vegandensity}
+\alias{densityplot.adonis}
+
+\title{
+  Kernel Density Estimation for Permutation Results in Vegan
+}
+
+\description{ 
+  The \code{density} functions can directly access the permutation
+  results of \pkg{vegan} functions, and \code{plot} can display the
+  densities. The \code{densityplot} method can access and display the
+  permutation results of functions that return permutations of several
+  statistics simultaneously.  
+}
+
+\usage{
+\method{density}{adonis}(x, ...)
+\method{plot}{vegandensity}(x, main = NULL, xlab = NULL, ylab = "Density", 
+   type = "l", zero.line = TRUE, obs.line = TRUE, ...)
+}
+
+\arguments{
+  \item{x}{The object to be handled. For \code{density} and
+     \code{densityplot} this is an object containing permutations. For
+     \code{plot} this is a result of \pkg{vegan} \code{density}
+     function.}
+  \item{main, xlab, ylab, type, zero.line}{Arguments of
+    \code{\link{plot.density}} and \code{\link[lattice]{densityplot}}
+    functions.}
+  \item{obs.line}{Draw vertical line for the observed
+    statistic. Logical value \code{TRUE} draws a red line, and
+    \code{FALSE} draws nothing. Alternatively, \code{obs.line} can be a
+    definition of the colour used for the line, either as a numerical
+    value from the \code{\link[grDevices]{palette}} or as the name of
+    the colour, or other normal definition of the colour.}
+  \item{\dots}{ Other arguments passed to the function. In
+    \code{density} these are passed to \code{\link{density.default}}.}
+}
+
+\details{ 
+
+  The \code{density} and \code{densityplot} function can directly access
+  permutation results of most \pkg{vegan} functions.  The \code{density}
+  function is identical to \code{\link{density.default}} and takes all
+  its arguments, but adds the observed statistic to the result as item
+  \code{"observed"}. The observed statistic is also put among the
+  permuted values so that the results are consistent with significance
+  tests. The \code{plot} method is similar to the default
+  \code{\link{plot.density}}, but can also add the observed statistic to
+  the graph as a vertical line.  The \code{densityplot} function is
+  based on the same function in the \pkg{lattice} package (see
+  \code{\link[lattice]{densityplot}}).
+
+  The density methods are available for \pkg{vegan} functions
+  \code{\link{adonis}}, \code{\link{anosim}}, \code{\link{mantel}},
+  \code{\link{mantel.partial}}, \code{\link{mrpp}},
+  \code{\link{permutest.cca}}, and \code{\link{protest}}.  The
+  \code{density} function for \code{\link{oecosimu}} is documented
+  separately, and it is also used for \code{\link{adipart}},
+  \code{\link{hiersimu}} and \code{\link{multipart}}.
+
+  All \pkg{vegan} \code{density} functions return an object of class
+  \code{"vegandensity"} inheriting from \code{\link{density}}, and can
+  be plotted with its \code{plot} method.  This is identical to the
+  standard \code{plot} of \code{densiy} objects, but can also add a
+  vertical line for the observed statistic.
+
+  Functions that can return several permuted statistics simultaenously
+  also have \code{\link[lattice]{densityplot}} method
+  (\code{\link{adonis}}, \code{\link{oecosimu}} and diversity partioning
+  functions based on \code{oecosimu}).  The standard
+  \code{\link{density}} can only handle univariate data, and a warning
+  is issued if the function is used for a model with several observed
+  statistics.  The \code{\link[lattice]{densityplot}} method is available
+  for \code{\link{adonis}} and \code{\link{oecosimu}} (documented
+  separately). NB, there is no \code{density} method for
+  \code{\link{anova.cca}}, but only for \code{\link{permutest.cca}}.
+
+}
+
+\value{
+  The \code{density} function returns the standard \code{\link{density}}
+  result object with one new item: \code{"observed"} for the observed
+  value of the statistic. The functions have a specific \code{plot}
+  method, but otherwise they use methods for
+  \code{\link{density.default}}, such as \code{print} and \code{lines}.
+}
+
+\author{
+  Jari Oksanen
+}
+
+\seealso{
+  \code{\link{density.default}}.
+}
+
+\examples{
+data(dune)
+data(dune.env)
+mod <- adonis(dune ~ Management, data = dune.env)
+plot(density(mod))
+library(lattice)
+mod <- adonis(dune ~ Management * Moisture, dune.env)
+densityplot(mod)
+}
+
+\keyword{ distribution }
+\keyword{ smooth }
diff --git a/man/designdist.Rd b/man/designdist.Rd
index b274c71..2f5734c 100644
--- a/man/designdist.Rd
+++ b/man/designdist.Rd
@@ -29,8 +29,8 @@ designdist(x, method = "(A+B-2*J)/(A+B)",
   \item{terms}{How shared and total components are found. For vectors
     \code{x} and \code{y} the  \code{"quadratic"} terms are \code{J = sum(x*y)},
     \code{A = sum(x^2)}, \code{B = sum(y^2)}, and \code{"minimum"} terms
-    are \code{J = sum(pmin(x,y))}, \code{A = sum(x)} and \code{B =
-      sum(y)}, and \code{"binary"} terms are either of these after transforming
+    are \code{J = sum(pmin(x,y))}, \code{A = sum(x)} and \code{B = sum(y)}, 
+    and \code{"binary"} terms are either of these after transforming
     data into binary form (shared number of species, and number of
     species for each row). }
   \item{abcd}{Use 2x2 contingency table notation for binary data:
@@ -65,15 +65,15 @@ designdist(x, method = "(A+B-2*J)/(A+B)",
   The function \code{designdist} can implement most dissimilarity
   indices in \code{\link{vegdist}} or elsewhere, and it can also be
   used to implement many other indices, amongst them, most of those
-  described in Legendre & Legendre (1998). It can also be used to
+  described in Legendre & Legendre (2012). It can also be used to
   implement all indices of beta diversity described in Koleff et
   al. (2003), but there also is a specific function
   \code{\link{betadiver}} for the purpose.
 
   If you want to implement binary dissimilarities based on the 2x2
   contingency table notation, you can set \code{abcd = TRUE}. In this
-  notation \code{a = J}, \code{b = A-J}, \code{c = B-J}, \code{d =
-  P-A-B+J}. This notation is often used instead of the more more
+  notation \code{a = J}, \code{b = A-J}, \code{c = B-J}, \code{d = P-A-B+J}. 
+  This notation is often used instead of the more more
   tangible default notation for reasons that are opaque to me. 
 }
 
@@ -85,7 +85,7 @@ designdist(x, method = "(A+B-2*J)/(A+B)",
   diversity for presence--absence data. \emph{J. Animal Ecol.}
   \strong{72}, 367--382. 
   
-  Legendre, P. and Legendre, L. (1998) \emph{Numerical Ecology}. 2nd
+  Legendre, P. and Legendre, L. (2012) \emph{Numerical Ecology}. 3rd
   English ed. Elsevier
   }
 \author{ Jari Oksanen }
diff --git a/man/envfit.Rd b/man/envfit.Rd
index da04ecc..c0768c2 100644
--- a/man/envfit.Rd
+++ b/man/envfit.Rd
@@ -6,6 +6,7 @@
 \alias{factorfit}
 \alias{plot.envfit}
 \alias{scores.envfit}
+\alias{labels.envfit}
 
 \title{Fits an Environmental Vector or Factor onto an Ordination }
 \description{
@@ -18,8 +19,8 @@
 \method{envfit}{default}(ord, env, permutations = 999, strata, choices=c(1,2), 
    display = "sites", w  = weights(ord), na.rm = FALSE, ...)
 \method{envfit}{formula}(formula, data, ...)
-\method{plot}{envfit}(x, choices = c(1,2), arrow.mul, at = c(0,0), axis = FALSE, 
-    p.max = NULL, col = "blue", bg, add = TRUE, ...)
+\method{plot}{envfit}(x, choices = c(1,2), labels, arrow.mul, at = c(0,0), 
+   axis = FALSE, p.max = NULL, col = "blue", bg, add = TRUE, ...)
 \method{scores}{envfit}(x, display, choices, ...)
 vectorfit(X, P, permutations = 0, strata, w, ...)
 factorfit(X, P, permutations = 0, strata, w, ...)
@@ -45,6 +46,13 @@ factorfit(X, P, permutations = 0, strata, w, ...)
     \code{na.rm = TRUE}.}
   \item{x}{A result object from \code{envfit}.}
   \item{choices}{Axes to plotted.}
+  \item{labels}{Change plotting labels. The argument should be a list
+    with elements \code{vectors} and \code{factors} which give the new
+    plotting labels. If either of these elements is omitted, the
+    default labels will be used. If there is only one type of elements
+    (only \code{vectors} or only \code{factors}), the labels can be
+    given as vector. The default labels can be displayed with
+    \code{labels} command.}
   \item{arrow.mul}{Multiplier for vector lengths. The arrows are
     automatically scaled similarly as in \code{\link{plot.cca}} if this
     is not given and \code{add = TRUE}.}
@@ -210,6 +218,11 @@ plot(ord, type = "n")
 ordispider(ord, Moisture, col="skyblue")
 points(ord, display = "sites", col = as.numeric(Moisture), pch=16)
 plot(fit, cex=1.2, axis=TRUE, bg = rgb(1, 1, 1, 0.5))
+## Use shorter labels for factor centroids
+labels(fit)
+plot(ord)
+plot(fit, labels=list(factors = paste("M", c(1,2,4,5), sep = "")),
+   bg = rgb(1,1,0,0.5))
 }
 \keyword{multivariate }
 \keyword{aplot}
diff --git a/man/fisherfit.Rd b/man/fisherfit.Rd
index 1df4768..ba462f4 100644
--- a/man/fisherfit.Rd
+++ b/man/fisherfit.Rd
@@ -111,8 +111,8 @@ as.preston(x, tiesplit = TRUE, ...)
   second octave are transferred to the higher one as well, but this is
   usually not as large a number of species. This practise makes data
   look more lognormal by reducing the usually high lowest
-  octaves. This can be achieved by setting argument \code{tiesplit =
-  TRUE}. With \code{tiesplit = FALSE} the frequencies are not split,
+  octaves. This can be achieved by setting argument \code{tiesplit = TRUE}. 
+  With \code{tiesplit = FALSE} the frequencies are not split,
   but all ones are in the lowest octave, all twos in the second, etc.
   Williamson & Gaston (2005) discuss alternative definitions in
   detail, and they should be consulted for a critical review of
@@ -128,8 +128,8 @@ as.preston(x, tiesplit = TRUE, ...)
   on the left so that some rare species are not observed. Function
   \code{prestonfit} fits the truncated lognormal model as a second
   degree log-polynomial to the octave pooled data using Poisson (when
-  \eqn{tiesplit = FALSE}) or quasi-Poisson (when \eqn{tiesplit =
-  TRUE}).  error. Function \code{prestondistr} fits left-truncated
+  \code{tiesplit = FALSE}) or quasi-Poisson (when \code{tiesplit = TRUE})
+  error.  Function \code{prestondistr} fits left-truncated
   Normal distribution to \eqn{\log_2}{log2} transformed non-pooled
   observations with direct maximization of log-likelihood. Function
   \code{prestondistr} is modelled after function
diff --git a/man/goodness.metaMDS.Rd b/man/goodness.metaMDS.Rd
index 65579d9..c2100b7 100644
--- a/man/goodness.metaMDS.Rd
+++ b/man/goodness.metaMDS.Rd
@@ -27,8 +27,8 @@
     number of points. }
   \item{p.col, l.col}{Point and line colours.}
   \item{lwd}{Line width. For \code{\link{monoMDS}} the default is
-    \code{lwd = 1} if more than two lines are drawn, and \code{lwd =
-    2} otherwise.}
+    \code{lwd = 1} if more than two lines are drawn, and \code{lwd = 2} 
+    otherwise.}
   \item{\dots}{Other parameters to functions, e.g. graphical parameters.}
 }
 \details{
diff --git a/man/humpfit.Rd b/man/humpfit.Rd
index adad855..9a9a75f 100644
--- a/man/humpfit.Rd
+++ b/man/humpfit.Rd
@@ -93,9 +93,9 @@ humpfit(mass, spno, family = poisson, start)
   methods. In addition, it can be accessed by the following methods
   for \code{\link{glm}} objects: \code{\link{AIC}},
   \code{\link{extractAIC}}, \code{\link{deviance}},
-  \code{\link{coef}}, \code{\link{residuals.glm}} (except \code{type =
-  "partial"}), \code{\link{fitted}}, and perhaps some others. In
-  addition, function \code{\link[ellipse]{ellipse.glm}} (package
+  \code{\link{coef}}, \code{\link{residuals.glm}} (except 
+  \code{type = "partial"}), \code{\link{fitted}}, and perhaps some others. 
+  In addition, function \code{\link[ellipse]{ellipse.glm}} (package
   \pkg{ellipse}) can be used to draw approximate confidence ellipses
   for pairs of parameters, if the normal assumptions look appropriate.
   } 
diff --git a/man/isomap.Rd b/man/isomap.Rd
index 755234a..96984e6 100644
--- a/man/isomap.Rd
+++ b/man/isomap.Rd
@@ -40,8 +40,8 @@ rgl.isomap(x, web = "white", ...)
 
   \item{type}{Plot observations either as \code{"points"},
     \code{"text"} or use \code{"none"} to plot no observations. The
-    \code{"text"} will use \code{\link{ordilabel}} if \code{net =
-    TRUE} and \code{\link{ordiplot}} if \code{net = FALSE}, and pass
+    \code{"text"} will use \code{\link{ordilabel}} if \code{net = TRUE} 
+    and \code{\link{ordiplot}} if \code{net = FALSE}, and pass
     extra arguments to these functions.}
 
   \item{web}{Colour of the web in \pkg{rgl} graphics.}
diff --git a/man/kendall.global.Rd b/man/kendall.global.Rd
index a7385cf..1d2d470 100644
--- a/man/kendall.global.Rd
+++ b/man/kendall.global.Rd
@@ -53,15 +53,16 @@ kendall.post(Y, group, nperm = 999, mult = "holm")
   compute Ward's agglomerative clustering of a matrix of correlations
   among the species. In detail: (1.1) compute a Pearson or Spearman
   correlation matrix (\code{correl.matrix}) among the species; (1.2)
-  turn it into a distance matrix: \code{mat.D =
-  as.dist(1-correl.matrix)}; (1.3) carry out Ward's hierarchical
-  clustering of that matrix using \code{hclust}: \code{clust.ward =
-  hclust(mat.D, "ward")}; (1.4) plot the dendrogram:
+  turn it into a distance matrix: \code{mat.D = as.dist(1-correl.matrix)}; 
+  (1.3) carry out Ward's hierarchical
+  clustering of that matrix using \code{hclust}: 
+  \code{clust.ward = hclust(mat.D, "ward")}; (1.4) plot the dendrogram:
   \code{plot(clust.ward, hang=-1)}; (1.5) cut the dendrogram in two
-  groups, retrieve the vector of species membership: \code{group.2 =
-  cutree(clust.ward, k=2)}. (1.6) After steps 2 and 3 below, you may
-  have to come back and try divisions of the species into k = 3, 4, 5,
-  \dots groups.
+  groups, retrieve the vector of species membership: 
+  \code{group.2 = cutree(clust.ward, k=2)}. (1.6) After steps 2 and 3 below, 
+  you may
+  have to come back and try divisions of the species into k = \eqn{3, 4, 5, \dots} 
+  groups.
 
   (2) Compute global tests of significance of the 2 (or more) groups
   using the function \code{kendall.global} and the vector defining the
diff --git a/man/make.cepnames.Rd b/man/make.cepnames.Rd
index 9637143..4ffca95 100644
--- a/man/make.cepnames.Rd
+++ b/man/make.cepnames.Rd
@@ -28,8 +28,8 @@ make.cepnames(names, seconditem = FALSE)
   function first makes valid \R names using \code{\link{make.names}},
   and then splits these into elemets. The CEP name is made by taking
   the four first letters of the first element, and four first letters
-  of the last (default) or the second element (with \code{seconditem =
-  TRUE}). If there was only one name element, it is
+  of the last (default) or the second element (with 
+  \code{seconditem = TRUE}). If there was only one name element, it is
   \code{\link{abbreviate}}d to eight letters. Finally, the names are
   made unique which may add numbers to duplicated names.
   
diff --git a/man/mantel.Rd b/man/mantel.Rd
index f5f109e..5357df9 100644
--- a/man/mantel.Rd
+++ b/man/mantel.Rd
@@ -13,9 +13,10 @@
 
 }
 \usage{
-mantel(xdis, ydis, method="pearson", permutations=999, strata)
+mantel(xdis, ydis, method="pearson", permutations=999, strata,
+    na.rm = FALSE)
 mantel.partial(xdis, ydis, zdis, method = "pearson", permutations = 999, 
-    strata)
+    strata, na.rm = FALSE)
 }
 
 \arguments{
@@ -26,6 +27,10 @@ mantel.partial(xdis, ydis, zdis, method = "pearson", permutations = 999,
   \item{strata}{An integer vector or factor specifying the strata for
     permutation. If supplied, observations are permuted only within the
     specified strata.}
+  \item{na.rm}{Remove missing values in calculation of Mantel
+    correlation. Use this option with care: Permutation tests can
+    be biased, in particular if two matrices had missing values in
+    matching positions.}
 }
 \details{
   Mantel statistic is simply a correlation between entries of two
@@ -56,21 +61,23 @@ mantel.partial(xdis, ydis, zdis, method = "pearson", permutations = 999,
     \code{\link{cor.test}}.}
   \item{statistic}{The Mantel statistic.}
   \item{signif}{Empirical significance level from permutations.}
-  \item{perm}{A vector of permuted values.}
+  \item{perm}{A vector of permuted values. The distribution of
+    permuted values can be inspected with \code{\link{density.mantel}} 
+    function.}
   \item{permutations}{Number of permutations.}
 }
 \references{ The test is due to Mantel, of course, but the
   current implementation is based on Legendre and Legendre.
 
-  Legendre, P. and Legendre, L. (1998) \emph{Numerical Ecology}. 2nd English
+  Legendre, P. and Legendre, L. (2012) \emph{Numerical Ecology}. 3rd English
   Edition. Elsevier.
   
 }
 
 \note{
-  Legendre & Legendre (1998) say that partial Mantel correlations 
-  often are difficult to interpret. 
-  }
+  Legendre & Legendre (2012, Box 10.4) warn against using partial
+  Mantel correlations.
+}
 
 \author{Jari Oksanen }
 
diff --git a/man/mantel.correlog.Rd b/man/mantel.correlog.Rd
index 2387ce8..745e4f4 100644
--- a/man/mantel.correlog.Rd
+++ b/man/mantel.correlog.Rd
@@ -7,8 +7,8 @@
 \description{
   Function \code{mantel.correlog} computes a multivariate
   Mantel correlogram. Proposed by Sokal (1986) and Oden and Sokal
-  (1986), the method is also described in Legendre and Legendre (1998,
-  pp. 736-738).
+  (1986), the method is also described in Legendre and Legendre (2012,
+  pp. 819--821).
 }
 
 \usage{
@@ -119,7 +119,7 @@ cutoff=TRUE, r.type="pearson", nperm=999, mult="holm", progressive=TRUE)
 
 \references{
 
-  Legendre, P. and L. Legendre. 1998. Numerical ecology, 2nd English
+  Legendre, P. and L. Legendre. 2012. Numerical ecology, 3rd English
   edition. Elsevier Science BV, Amsterdam.
 
   Mantel, N. 1967. The detection of disease clustering and a generalized
@@ -146,10 +146,10 @@ mite.hel <- decostand(mite, "hellinger")
 mite.hel.resid <- resid(lm(as.matrix(mite.hel) ~ ., data=mite.xy))
 
 # Compute the detrended species distance matrix
-mite.hel.D = dist(mite.hel.resid)
+mite.hel.D <- dist(mite.hel.resid)
 
 # Compute Mantel correlogram with cutoff, Pearson statistic
-mite.correlog = mantel.correlog(mite.hel.D, XY=mite.xy, nperm=49)
+mite.correlog <- mantel.correlog(mite.hel.D, XY=mite.xy, nperm=49)
 summary(mite.correlog)
 mite.correlog   
 # or: print(mite.correlog)
@@ -157,8 +157,8 @@ mite.correlog
 plot(mite.correlog)
 
 # Compute Mantel correlogram without cutoff, Spearman statistic
-mite.correlog2 = mantel.correlog(mite.hel.D, XY=mite.xy, cutoff=FALSE, 
-r.type="spearman", nperm=49)
+mite.correlog2 <- mantel.correlog(mite.hel.D, XY=mite.xy, cutoff=FALSE, 
+   r.type="spearman", nperm=49)
 summary(mite.correlog2)
 mite.correlog2
 plot(mite.correlog2)
diff --git a/man/metaMDS.Rd b/man/metaMDS.Rd
index ecddf2e..93c26a8 100644
--- a/man/metaMDS.Rd
+++ b/man/metaMDS.Rd
@@ -114,8 +114,7 @@ metaMDSredist(object, ...)
    \code{"fail"} or \code{"add"} a small positive value, or
    \code{"ignore"}. \code{\link{monoMDS}} accepts zero dissimilarities
    and the default is \code{zerodist = "ignore"}, but with
-   \code{\link[MASS]{isoMDS}} you may need to set \code{zerodist =
-   "add"}.}
+   \code{\link[MASS]{isoMDS}} you may need to set \code{zerodist = "add"}.}
 
  \item{distfun}{Dissimilarity function. Any function returning a
    \code{dist} object and accepting argument \code{method} can be used
diff --git a/man/mrpp.Rd b/man/mrpp.Rd
index 85344a4..b874e40 100644
--- a/man/mrpp.Rd
+++ b/man/mrpp.Rd
@@ -138,7 +138,8 @@ The function returns a list of class mrpp with following items:
     the dist object.}
   \item{weight.type}{The choice of group weights used.}
   \item{boot.deltas}{The vector of "permuted deltas," the deltas
-    calculated from each of the permuted datasets.}
+    calculated from each of the permuted datasets. The distribution of
+    this item can be inspected with \code{\link{density.mrpp}} function.}
   \item{permutations}{The number of permutations used.}
 }
 \references{
diff --git a/man/multipart.Rd b/man/multipart.Rd
index aaf540a..eace009 100644
--- a/man/multipart.Rd
+++ b/man/multipart.Rd
@@ -3,7 +3,6 @@
 \alias{multipart}
 \alias{multipart.default}
 \alias{multipart.formula}
-\alias{print.multipart}
 \title{Multiplicative Diversity Partitioning}
 
 \description{
@@ -26,10 +25,9 @@ multipart(...)
     all rows are in the same group in the second level.}
   \item{formula}{A two sided model formula in the form \code{y ~ x}, where \code{y} 
     is the community data matrix with samples as rows and species as column. Right 
-    hand side (\code{x}) must contain factors referring to levels of sampling hierarchy, 
+    hand side (\code{x}) must grouping vaiables referring to levels of sampling hierarchy, 
     terms from right to left will be treated as nested (first column is the lowest, 
-    last is the highest level). These variables must be factors in order to unambiguous 
-    handling. Interaction terms are not allowed.}
+    last is the highest level, at least two levels specified). Interaction terms are not allowed.}
   \item{data}{A data frame where to look for variables defined in the right hand side 
     of \code{formula}. If missing, variables are looked in the global environment.}
   \item{index}{Character, the entropy index to be calculated (see Details).}
@@ -105,10 +103,10 @@ cutter <- function (x, cut = seq(0, 10, by = 2.5)) {
     out <- rep(1, length(x))
     for (i in 2:(length(cut) - 1))
         out[which(x > cut[i] & x <= cut[(i + 1)])] <- i
-    return(as.factor(out))}
+    return(out)}
 ## The hierarchy of sample aggregation
 levsm <- data.frame(
-    l1=as.factor(1:nrow(mite)),
+    l1=1:nrow(mite),
     l2=cutter(mite.xy$y, cut = seq(0, 10, by = 2.5)),
     l3=cutter(mite.xy$y, cut = seq(0, 10, by = 5)),
     l4=cutter(mite.xy$y, cut = seq(0, 10, by = 10)))
diff --git a/man/ordiplot3d.Rd b/man/ordiplot3d.Rd
index 64c956b..fccb920 100644
--- a/man/ordiplot3d.Rd
+++ b/man/ordiplot3d.Rd
@@ -41,8 +41,7 @@ orglspider(object, groups, display = "sites", w = weights(object, display),
     text labels.}
   \item{ax.col}{Axis colour (concerns only the crossed axes through the
     origin).}
-  \item{text}{Text to override the default with \code{type =
-    "t"}.}
+  \item{text}{Text to override the default with \code{type = "t"}.}
   \item{envfit}{Fitted environmental variables from \code{\link{envfit}}
     displayed in the graph.}
   \item{xlab, ylab, zlab}{Axis labels passed to
diff --git a/man/ordisurf.Rd b/man/ordisurf.Rd
index 1c079a9..2260d4c 100644
--- a/man/ordisurf.Rd
+++ b/man/ordisurf.Rd
@@ -113,9 +113,9 @@
   Wood, 2011). The addition of this extra penalty is invoked by
   setting argument \code{select} to \code{TRUE}. The function plots
   the fitted contours with convex hull of data points either over an
-  existing ordination diagram or draws a new plot. If \code{select ==
-  TRUE} and the smooth is effectively penalised out of the model, no
-  contours will be plotted.
+  existing ordination diagram or draws a new plot. If 
+  \code{select = TRUE} and the smooth is effectively penalised out of 
+  the model, no contours will be plotted.
 
   \code{\link[mgcv]{gam}} determines the degree of smoothness for the
   fitted response surface during model fitting. Argument \code{method}
diff --git a/man/predict.cca.Rd b/man/predict.cca.Rd
index cf124ad..b522497 100644
--- a/man/predict.cca.Rd
+++ b/man/predict.cca.Rd
@@ -44,8 +44,8 @@
   
   \item{newdata}{New data frame to be used in prediction or in
     calibration.  Usually this a new community data frame, but with
-    \code{type = "lc"} and for constrained component with \code{type =
-    "response"} and \code{type = "working"} it must be an environment
+    \code{type = "lc"} and for constrained component with 
+    \code{type = "response"} and \code{type = "working"} it must be a 
     data frame.  The \code{newdata} must have the same number of rows as
     the original community data for a \code{\link{cca}} result with
     \code{type = "response"} or \code{type = "working"}.  If the
@@ -112,9 +112,9 @@
   \code{type = "lc"} the function finds the linear combination scores
   for sites from environmental data. In that case the new data frame
   must contain all constraining and conditioning environmental variables
-  of the model formula. With \code{type = "response"} or \code{type =
-  "working"} the new data must contain envinronmental variables if
-  constrained component is desired, and community data matrix if
+  of the model formula. With \code{type = "response"} or 
+  \code{type = "working"} the new data must contain envinronmental variables 
+  if constrained component is desired, and community data matrix if
   residual or unconstrained component is desired.  With these types, the
   function uses \code{newdata} to find new \code{"lc"} (constrained) or
   \code{"wa"} scores (unconstrained) and then finding the response or
diff --git a/man/procrustes.Rd b/man/procrustes.Rd
index 03abe52..9d0a19c 100644
--- a/man/procrustes.Rd
+++ b/man/procrustes.Rd
@@ -9,7 +9,6 @@
 \alias{fitted.procrustes}
 \alias{predict.procrustes}
 \alias{protest}
-\alias{print.protest}
 
 \title{Procrustes Rotation of Two Configurations and PROTEST }
 \description{
@@ -161,7 +160,9 @@ protest(X, Y, scores = "sites", permutations = 999, strata, ...)
   \item{call}{Function call.}
   \item{t0}{This and the following items are only in class
     \code{protest}:  Procrustes correlation from non-permuted solution.}
-  \item{t}{Procrustes correlations from permutations.}
+  \item{t}{Procrustes correlations from permutations. The distribution
+    of these correlations can be inspected with \code{\link{density.protest}} 
+    function.}
   \item{signif}{`Significance' of \code{t}}
   \item{permutations}{Number of permutations.}
   \item{strata}{The name of the stratifying variable.}
diff --git a/man/radfit.Rd b/man/radfit.Rd
index 77427a2..8993334 100644
--- a/man/radfit.Rd
+++ b/man/radfit.Rd
@@ -3,22 +3,34 @@
 \alias{radfit.default}
 \alias{radfit.data.frame}
 \alias{AIC.radfit}
+\alias{AIC.radfit.frame}
 \alias{as.rad}
 \alias{coef.radfit}
+\alias{coef.radfit.frame}
+\alias{deviance.radfit}
+\alias{deviance.radfit.frame}
+\alias{logLik, radfit}
+\alias{logLik, radfit.frame}
 \alias{fitted.radfit}
+\alias{fitted.radfit.frame}
 \alias{lines.radline}
+\alias{lines.radfit}
 \alias{plot.radfit.frame}
 \alias{plot.radfit}
 \alias{plot.radline}
 \alias{plot.rad}
 \alias{radlattice}
 \alias{points.radline}
+\alias{points.radfit}
 \alias{summary.radfit.frame}
 \alias{rad.preempt}
 \alias{rad.lognormal}
 \alias{rad.zipf}
 \alias{rad.zipfbrot}
 \alias{rad.null}
+\alias{predict.radline}
+\alias{predict.radfit}
+\alias{predict.radfit.frame}
 
 \title{ Rank -- Abundance or Dominance / Diversity Models}
 \description{
@@ -27,152 +39,167 @@
   Zipf and Zipf-Mandelbrot models of species abundance.
 }
 \usage{
-\method{radfit}{data.frame}(df, ...)
-\method{plot}{radfit.frame}(x, order.by, BIC = FALSE, model, legend = TRUE,
-     as.table = TRUE, ...)
 \method{radfit}{default}(x, ...)
-\method{plot}{radfit}(x, BIC = FALSE, legend = TRUE, ...)  
-radlattice(x, BIC = FALSE, ...)
 rad.null(x, family=poisson, ...)
 rad.preempt(x, family = poisson, ...)
 rad.lognormal(x, family = poisson, ...)
 rad.zipf(x, family = poisson, ...)
 rad.zipfbrot(x, family = poisson, ...)
+\method{predict}{radline}(object, newdata, total, ...)
+\method{plot}{radfit}(x, BIC = FALSE, legend = TRUE, ...)  
+\method{plot}{radfit.frame}(x, order.by, BIC = FALSE, model, legend = TRUE,
+     as.table = TRUE, ...)
 \method{plot}{radline}(x, xlab = "Rank", ylab = "Abundance", type = "b", ...)
-\method{lines}{radline}(x, ...)
-\method{points}{radline}(x, ...)
+radlattice(x, BIC = FALSE, ...)
+\method{lines}{radfit}(x, ...)
+\method{points}{radfit}(x, ...)
 as.rad(x)
 \method{plot}{rad}(x, xlab = "Rank", ylab = "Abundance", log = "y", ...)
 }
 
 \arguments{
-  \item{df}{Data frame where sites are rows and species are columns.}
-  \item{x}{A vector giving species abundances in a site, or an object to
+  \item{x}{Data frame, matrix or a vector giving species abundances, or an object to
     be plotted.}
+
+  \item{family}{Error distribution (passed to \code{\link{glm}}). All
+    alternatives accepting \code{link = "log"} in \code{\link{family}}
+    can be used, although not all make sense.}
+
+  \item{object}{A fitted result object.}
+
+  \item{newdata}{Ranks used for ordinations. All models can
+    interpolate to non-integer \dQuote{ranks} (although this may be
+    approximate), but extrapolation may fail}
+
+  \item{total}{The new total used for predicting abundance. Observed
+    total count is used if this is omitted.}
+
   \item{order.by}{A vector used for ordering sites in plots.}
   \item{BIC}{Use Bayesian Information Criterion, BIC, instead of
-    Akaike's AIC. The penalty for a parameter is \eqn{k = \log(S)}{k =
+    Akaike's AIC. The penalty in BIC is \eqn{k = \log(S)}{k =
       log(S)}  where \eqn{S} is the number of species, whereas AIC uses
     \eqn{k = 2}.} 
-  \item{model}{Show only the specified model. If missing, AIC is used to
-    select the model. The model names (which can be abbreviated) are
-    \code{Preemption}, \code{Lognormal}, \code{Veiled.LN},
-    \code{Zipf}, \code{Mandelbrot}. }
+  \item{model}{Show only the specified model. If missing, AIC is used
+    to select the model. The model names (which can be abbreviated)
+    are \code{Null}, \code{Preemption}, \code{Lognormal}, \code{Zipf},
+    \code{Mandelbrot}. }
   \item{legend}{Add legend of line colours.}
   \item{as.table}{Arrange panels starting from upper left corner (passed
     to \code{\link[lattice]{xyplot}}).}
-  \item{family}{Error distribution (passed to \code{\link{glm}}). All
-    alternatives accepting \code{link = "log"} in \code{\link{family}}
-    can be used, although not all make sense.}
   \item{xlab,ylab}{Labels for \code{x} and \code{y} axes.}
   \item{type}{Type of the plot, \code{"b"} for plotting both observed points
     and fitted lines, \code{"p"} for only points, \code{"l"} for only
     fitted lines, and \code{"n"} for only setting the frame. }
   \item{log}{Use logarithmic scale for given axis. The default
-    \code{log =" y"} gives the traditional plot in community ecology
+    \code{log = "y"} gives the traditional plot of community ecology
     where the pre-emption model is a straight line, and with
     \code{log = "xy"} Zipf model is a straight line. With
     \code{log = ""} both axes are in the original arithmetic scale.}
   \item{\dots}{Other parameters to functions. }
 }
 \details{
-  Rank -- Abundance Dominance (RAD) or Dominance/Diversity plots
+
+  Rank--Abundance Dominance (RAD) or Dominance/Diversity plots
   (Whittaker 1965) display logarithmic species abundances against
-  species rank order. These plots are supposed to be
-  effective in analysing types of abundance distributions in
-  communities. These functions fit some of the most popular models mainly
-  following Wilson (1991). Function \code{as.rad} constructs observed
-  RAD data.
-  Functions \code{rad.XXXX} (where \code{XXXX} is a name) fit
-  the individual models, and
-  function \code{radfit} fits all models.  The
-  argument of the function \code{radfit} can be either a vector for a
-  single community or a data frame where each row represents a
-  distinct community.  All these functions have their own \code{plot}
-  functions. When the argument is a data frame, \code{plot} uses
-  \code{\link[lattice]{Lattice}} graphics, and other \code{plot} functions use
-  ordinary graphics. The ordinary graphics functions return invisibly an
-  \code{\link{ordiplot}} object for observed points, and function
-  \code{\link{identify.ordiplot}} can be used to label selected
-  species. The most complete control of graphics can be achieved
-  with \code{rad.XXXX} methods which have \code{points} and \code{lines}
-  functions to add observed values and fitted models into existing
-  graphs.  Alternatively, \code{radlattice} uses
-  \code{\link[lattice]{Lattice}} graphics to display each
-  \code{radfit} model in a separate panel together with their AIC or
-  BIC values.
+  species rank order. These plots are supposed to be effective in
+  analysing types of abundance distributions in communities. These
+  functions fit some of the most popular models mainly following
+  Wilson (1991).
 
+  Functions \code{rad.null}, \code{rad.preempt}, \code{rad.lognormal},
+  \code{rad.zipf} and \code{zipfbrot} fit the individual models
+  (described below) for a single vector (row of data frame), and
+  function \code{radfit} fits all models. The argument of the function
+  \code{radfit} can be either a vector for a single community or a data
+  frame where each row represents a distinct community.
+  
   Function \code{rad.null} fits a brokenstick model where the expected
-  abundance of species at rank \eqn{r} is \eqn{a_r = (J/S) \sum_{x=r}^S
-    (1/x)}{a[r] = J/S sum(from x=r to S) 1/x} (Pielou 1975), where \eqn{J}
-  is the total number of individuals (site total) and \eqn{S} is the
-  total number of species in the community.  This gives a Null model
-  where the individuals are randomly distributed among observed species,
-  and there are no fitted parameters. 
+  abundance of species at rank \eqn{r} is \eqn{a_r = (J/S)
+  \sum_{x=r}^S (1/x)}{a[r] = J/S sum(from x=r to S) 1/x} (Pielou
+  1975), where \eqn{J} is the total number of individuals (site total)
+  and \eqn{S} is the total number of species in the community.  This
+  gives a Null model where the individuals are randomly distributed
+  among observed species, and there are no fitted parameters.
   Function \code{rad.preempt} fits the niche preemption model,
   a.k.a. geometric series or Motomura model, where the expected
-  abundance \eqn{a} of species at rank \eqn{r} is \eqn{a_r = J \alpha (1 -
-    \alpha)^{r-1}}{a[r] = J*alpha*(1-alpha)^(r-1)}. The only estimated
-  parameter is the preemption coefficient \eqn{\alpha} which gives the
-  decay rate of abundance per rank.
-  The niche preemption model is a straight line in a
-  RAD plot. Function \code{rad.lognormal} fits a log-Normal model which
-  assumes that the logarithmic abundances are distributed Normally, or
-  \eqn{a_r =  \exp( \log \mu + \log \sigma N)}{a[r] = exp(log(mu) +
-    log(sigma) * N)}, where \eqn{N} is a Normal deviate. 
-  Function \code{rad.zipf} fits the Zipf model \eqn{a_r = J p_1
-    r^\gamma}{a[r] = J*p1*r^gamma} where \eqn{p_1}{p1} is the fitted
-  proportion of the most abundant species, and \eqn{\gamma} is a
-  decay coefficient. The
-  Zipf -- Mandelbrot 
-  model (\code{rad.zipfbrot}) adds one parameter: \eqn{a_r = J c
-    (r + \beta)^\gamma}{a[r] = J*c*(r+beta)^gamma} after which
-  \eqn{p_1}{p1} of the Zipf model changes into a meaningless scaling
-  constant \eqn{c}. There are grand narratives about ecological
-  mechanisms behind each model (Wilson 1991), but
-  several alternative and contrasting mechanisms can produce
-  similar models and a good fit does not imply a specific mechanism.
-
-  Log-Normal and Zipf models are generalized linear
-  models (\code{\link{glm}}) with logarithmic link function.
-  Zipf-Mandelbrot adds one
-  nonlinear parameter to the Zipf model, and is fitted using
+  abundance \eqn{a} of species at rank \eqn{r} is \eqn{a_r = J \alpha
+  (1 - \alpha)^{r-1}}{a[r] = J*alpha*(1-alpha)^(r-1)}. The only
+  estimated parameter is the preemption coefficient \eqn{\alpha} which
+  gives the decay rate of abundance per rank.  The niche preemption
+  model is a straight line in a RAD plot.  Function
+  \code{rad.lognormal} fits a log-Normal model which assumes that the
+  logarithmic abundances are distributed Normally, or \eqn{a_r = \exp(
+  \log \mu + \log \sigma N)}{a[r] = exp(log(mu) + log(sigma) * N)},
+  where \eqn{N} is a Normal deviate.  Function \code{rad.zipf} fits
+  the Zipf model \eqn{a_r = J p_1 r^\gamma}{a[r] = J*p1*r^gamma} where
+  \eqn{p_1}{p1} is the fitted proportion of the most abundant species,
+  and \eqn{\gamma} is a decay coefficient. The Zipf--Mandelbrot model
+  (\code{rad.zipfbrot}) adds one parameter: \eqn{a_r = J c (r +
+  \beta)^\gamma}{a[r] = J*c*(r+beta)^gamma} after which \eqn{p_1}{p1}
+  of the Zipf model changes into a meaningless scaling constant
+  \eqn{c}. 
+
+  Log-Normal and Zipf models are generalized linear models
+  (\code{\link{glm}}) with logarithmic link function.  Zipf--Mandelbrot
+  adds one nonlinear parameter to the Zipf model, and is fitted using
   \code{\link{nlm}} for the nonlinear parameter and estimating other
-  parameters and log-Likelihood with \code{\link{glm}}. Pre-emption
-  model is fitted as purely nonlinear model. There are no estimated
-  parameters in the Null model.  The default
-  \code{\link{family}} is \code{poisson} which is appropriate only for
-  genuine counts (integers), but other families that accept \code{link =
-    "log"} can be used. Family \code{\link{Gamma}} may be
-  appropriate for abundance data, such as cover. The ``best''
-  model is selected by \code{\link{AIC}}. Therefore ``quasi'' families
-  such as \code{\link{quasipoisson}} cannot be used: they do not
-  have \code{\link{AIC}} nor log-Likelihood needed in non-linear
-  models.
+  parameters and log-Likelihood with \code{\link{glm}}. Preemption
+  model is fitted as a purely nonlinear model. There are no estimated
+  parameters in the Null model.  
+
+  The default \code{\link{family}} is \code{poisson} which is
+  appropriate only for genuine counts (integers), but other families
+  that accept \code{link = "log"} can be used. Families
+  \code{\link{Gamma}} or \code{\link{gaussian}} may be appropriate for
+  abundance data, such as cover. The ``best'' model is selected by
+  \code{\link{AIC}}. Therefore ``quasi'' families such as
+  \code{\link{quasipoisson}} cannot be used: they do not have
+  \code{\link{AIC}} nor log-Likelihood needed in non-linear models.
+  
+  All these functions have their own \code{plot} functions. When
+  \code{radfit} was applied for a data frame, \code{plot} uses
+  \code{\link[lattice]{Lattice}} graphics, and other \code{plot}
+  functions use ordinary graphics. The ordinary graphics functions
+  return invisibly an \code{\link{ordiplot}} object for observed points,
+  and function \code{\link{identify.ordiplot}} can be used to label
+  selected species.  Alternatively, \code{radlattice} uses
+  \code{\link[lattice]{Lattice}} graphics to display each \code{radfit}
+  model of a single site in a separate panel together with their AIC or
+  BIC values.
+
+  Function \code{as.rad} is a base function to construct ordered RAD
+  data. Its \code{plot} is used by other RAD \code{plot} functions
+  which pass extra arguments (such as \code{xlab} and \code{log}) to
+  this function.
+
 }
 
 \value{
-  Function \code{rad.XXXX} will return an object of class
-  \code{radline}, which is constructed to resemble results of \code{\link{glm}}
-  and has many (but not all) of its components, even when only
-  \code{\link{nlm}} was used in fitting. At least the following
-  \code{\link{glm}} methods can be applied to the result:
-  \code{\link{fitted}}, \code{\link{residuals.glm}}  with alternatives
-  \code{"deviance"} (default), \code{"pearson"}, \code{"response"},
-  function \code{\link{coef}}, \code{\link{AIC}},
-  \code{\link{extractAIC}}, and \code{\link{deviance}}.
-  Function \code{radfit} applied to a vector will return
-  an object of class \code{radfit} with item \code{y} for the
-  constructed RAD, item \code{family} for the error distribution, and
-  item \code{models} containing each \code{radline} object as an
-  item. In addition, there are special \code{AIC}, \code{coef} and
-  \code{fitted} implementations for \code{radfit} results. 
-  When applied to a data frame
-  \code{radfit} will return an object of class \code{radfit.frame} which
-  is a list of \code{radfit} objects; function \code{summary} can be
-  used to display the results for individual \code{radfit} objects.
-  The functions are still
-  preliminary, and the items in the \code{radline} objects may change.
+
+  Functions \code{rad.null}, \code{rad.preempt}, \code{rad.lognormal},
+  \code{zipf} and \code{zipfbrot} fit each a single RAD model to a
+  single site. The result object has class \code{"radline"} and
+  inherits from \code{\link{glm}}, and can be handled by some (but not
+  all) \code{\link{glm}} methods.
+
+  Function \code{radfit} fits all models either to a single site or to
+  all rows of a data frame or a matrix. When fitted to a single site,
+  the function returns an object of class \code{"radfit"} with items
+  \code{y} (observed values), \code{\link{family}}, and \code{models}
+  which is a list of fitted \code{"radline"} models.  When applied for a
+  data frame or matrix, \code{radfit} function returns an object of
+  class \code{"radfit.frame"} which is a list of \code{"radfit"}
+  objects, each item names by the corresponding row name.
+
+  All result objects (\code{"radline"}, \code{"radfit"},
+  \code{"radfit.frame"}) can be accessed with same method functions.
+  The following methods are available: \code{\link{AIC}},
+  \code{\link{coef}}, \code{\link{deviance}}, \code{\link{logLik}}. In
+  addition the fit results can be accessed with \code{\link{fitted}},
+  \code{\link{predict}} and \code{\link{residuals}} (inheriting from
+  \code{\link{residuals.glm}}).The graphical functions were discussed
+  above in Details.
+
 }
 
 \references{
diff --git a/man/raupcrick.Rd b/man/raupcrick.Rd
index 05d0aeb..5767755 100644
--- a/man/raupcrick.Rd
+++ b/man/raupcrick.Rd
@@ -85,8 +85,8 @@ raupcrick(comm, null = "r1", nsimul = 999, chase = FALSE)
   \pkg{vegan}.  Chase et al. (2011) script with \code{split = TRUE}
   uses half of tied simulation values to calculate a distance measure,
   and that choice cannot be directly reproduced in vegan (it is the
-  average of \pkg{vegan} \code{raupcrick} results with \code{chase =
-  TRUE} and \code{chase = FALSE}).}
+  average of \pkg{vegan} \code{raupcrick} results with 
+  \code{chase = TRUE} and \code{chase = FALSE}).}
 
 \seealso{The function is based on \code{\link{oecosimu}}. Function
   \code{\link{vegdist}} with {method = "raup"} implements a related
diff --git a/man/screeplot.cca.Rd b/man/screeplot.cca.Rd
index 2e06711..ee83d07 100644
--- a/man/screeplot.cca.Rd
+++ b/man/screeplot.cca.Rd
@@ -81,7 +81,7 @@ bstick(n, \dots)
 
   Function \code{bstick} gives the brokenstick values which are ordered
   random proportions, defined as  \eqn{p_i = (tot/n) \sum_{x=i}^n 
-    (1/x)}{p[i] = tot/n sum(from x=i to n) 1/x} (Legendre & Legendre 1998), where
+    (1/x)}{p[i] = tot/n sum(from x=i to n) 1/x} (Legendre & Legendre 2012), where
   \eqn{tot} is the total  and \eqn{n} is the number of brokenstick
   components (cf. \code{\link{radfit}}).  Broken stick has
   been recommended as a stopping rule in principal component analysis
@@ -113,7 +113,7 @@ bstick(n, \dots)
   analysis: a comparison of heuristical and statistical
   approaches. \emph{Ecology} 74, 2204--2214.
 
-  Legendre, P. and Legendre, L. (1998) \emph{Numerical Ecology}. 2nd English
+  Legendre, P. and Legendre, L. (2012) \emph{Numerical Ecology}. 3rd English
   ed. Elsevier.
   }
 \note{Function \code{screeplot} is generic from \code{R} version
diff --git a/man/specaccum.Rd b/man/specaccum.Rd
index 6ecff3a..9beacf6 100644
--- a/man/specaccum.Rd
+++ b/man/specaccum.Rd
@@ -39,8 +39,7 @@ fitspecaccum(object, model, method = "random", ...)
     expected richness following
     Coleman et al. 1982, and \code{"rarefaction"} finds the mean when
     accumulating individuals instead of sites.  }
-  \item{permutations}{Number of permutations with \code{method =
-      "random"}.}
+  \item{permutations}{Number of permutations with \code{method = "random"}.}
   \item{conditioned}{ Estimation of standard deviation is conditional on
     the empirical dataset for the exact SAC}
   \item{gamma}{Method for estimating the total extrapolated number of species in the
@@ -101,8 +100,8 @@ fitspecaccum(object, model, method = "random", ...)
   achieves this by applying function \code{\link{rarefy}} with number of individuals
   corresponding to average number of individuals per site.
 
-  The function has a \code{plot} method. In addition, \code{method =
-  "random"} has \code{summary} and \code{boxplot} methods. 
+  The function has a \code{plot} method. In addition, \code{method = "random"} 
+  has \code{summary} and \code{boxplot} methods. 
 
   Function \code{predict} can return the values corresponding to
   \code{newdata} using linear (\code{\link{approx}}) or spline
diff --git a/man/stepacross.Rd b/man/stepacross.Rd
index c3ba571..0462a46 100644
--- a/man/stepacross.Rd
+++ b/man/stepacross.Rd
@@ -46,8 +46,8 @@ stepacross(dis, path = "shortest", toolong = 1, trace = TRUE, ...)
   \code{toolong}. 
 
   De'ath (1999) suggested a simplified method known as extended
-  dissimilarities, which are calculated with \code{path =
-    "extended"}. In this method, dissimilarities that are
+  dissimilarities, which are calculated with \code{path = "extended"}. 
+  In this method, dissimilarities that are
   \code{toolong} or longer are first made \code{NA}, and then the function
   tries to replace these \code{NA} dissimilarities with a path through
   single stepping stone points. If not all \code{NA} could be 
@@ -68,8 +68,8 @@ stepacross(dis, path = "shortest", toolong = 1, trace = TRUE, ...)
   \code{toolong} is zero or negative, the function does not make any
   dissimilarities into \code{NA}. If there are no \code{NA}s in the
   input  and \code{toolong = 0}, \code{path = "shortest"}
-  will find shorter paths for semi-metric indices, and \code{path =
-    "extended"} will do nothing. Function \code{\link{no.shared}} can be
+  will find shorter paths for semi-metric indices, and \code{path = "extended"} 
+  will do nothing. Function \code{\link{no.shared}} can be
   used to set dissimilarities to \code{NA}.
   
   If the data are disconnected or there is no path between all points,
@@ -82,8 +82,8 @@ stepacross(dis, path = "shortest", toolong = 1, trace = TRUE, ...)
 
   Alternative \code{path = "shortest"} uses Dijkstra's method for
   finding flexible shortest paths, implemented as priority-first search
-  for dense graphs (Sedgewick 1990). Alternative \code{path =
-    "extended"} follows De'ath (1999), but implementation is simpler
+  for dense graphs (Sedgewick 1990). Alternative \code{path = "extended"} 
+  follows De'ath (1999), but implementation is simpler
   than in his code.
   
 }
diff --git a/man/treedive.Rd b/man/treedive.Rd
index 24a0cf7..82e7d56 100644
--- a/man/treedive.Rd
+++ b/man/treedive.Rd
@@ -66,9 +66,9 @@ treedist(x, tree, relative = TRUE, ...)
   species traits contain \code{\link{factor}} or \code{\link{ordered}}
   factor variables, it is recommended to use Gower distances for mixed
   data (function \code{\link[cluster]{daisy}} in package \pkg{cluster}),
-  and usually the recommended clustering method is UPGMA (\code{method =
-  "average"} in function \code{\link{hclust}}) (Podani and Schmera
-  2006).
+  and usually the recommended clustering method is UPGMA 
+  (\code{method = "average"} in function \code{\link{hclust}}) 
+  (Podani and Schmera 2006).
 
   It is possible to analyse the non-randomness of functional diversity
   using \code{\link{oecosimu}}. This needs specifying an adequate Null
diff --git a/man/varpart.Rd b/man/varpart.Rd
index babe89a..edee8e9 100644
--- a/man/varpart.Rd
+++ b/man/varpart.Rd
@@ -168,7 +168,7 @@ table. }
 Borcard, D., P. Legendre & P. Drapeau. 1992. Partialling out the spatial
 component of ecological variation. Ecology 73: 1045--1055.
 
-Legendre, P. & L. Legendre. 1998. Numerical ecology, 2nd English edition.
+Legendre, P. & L. Legendre. 2012. Numerical ecology, 3rd English edition.
 Elsevier Science BV, Amsterdam.
 
 (b) Reference on transformations for species data
@@ -228,8 +228,12 @@ vegandocs("partition")
 # Formula shortcut "~ ." means: use all variables in 'data'.
 mod <- varpart(mite, ~ ., mite.pcnm, data=mite.env, transfo="hel")
 mod
-showvarparts(2)
-plot(mod)
+
+## argument 'bg' is passed to circle drawing, and the following
+## defines semitransparent colours
+col2 <- rgb(c(1,1),c(1,0), c(0,1), 0.3)
+showvarparts(2, bg = col2)
+plot(mod, bg = col2)
 # Alternative way of to conduct this partitioning
 # Change the data frame with factors into numeric model matrix
 mm <- model.matrix(~ SubsDens + WatrCont + Substrate + Shrub + Topo, mite.env)[,-1]
diff --git a/man/vegan-internal.Rd b/man/vegan-internal.Rd
index e28e105..1bd9666 100644
--- a/man/vegan-internal.Rd
+++ b/man/vegan-internal.Rd
@@ -10,6 +10,7 @@
 \alias{ordiArrowMul}
 \alias{ordiArgAbsorber}
 \alias{veganCovEllipse}
+\alias{hierParseFormula}
 
 \title{Internal vegan functions}
 
@@ -31,6 +32,7 @@ centroids.cca(x, mf, wt)
 permuted.index(n, strata)
 pasteCall(call, prefix = "Call:")
 veganCovEllipse(cov, center = c(0, 0), scale = 1, npoints = 100)
+hierParseFormula(formula, data)
 }
 
 \details{ The description here is only intended for \pkg{vegan}
@@ -79,6 +81,11 @@ veganCovEllipse(cov, center = c(0, 0), scale = 1, npoints = 100)
 
   \code{veganCovEllipse} finds the coordinates for drawing a
   covariance ellipse.
+
+  \code{hierParseFormula} returns a list of one matrix (left hand side)
+  and a model frame with factors representing hierarchy levels 
+  (right hand side) to be used in \code{\link{adipart}}, 
+  \code{\link{multipart}} and \code{\link{hiersimu}}.
 }
 
 \keyword{internal }
diff --git a/man/vegdist.Rd b/man/vegdist.Rd
index 2206ed8..2a19cb3 100644
--- a/man/vegdist.Rd
+++ b/man/vegdist.Rd
@@ -170,26 +170,24 @@
   are divided by \eqn{\log(2)}{log(2)} so that the results will be in
   the conventional range \eqn{0 \dots 1}.
 
-  Raup--Crick dissimilarity (\code{method = "raup"}) is a
-  probabilistic index based on presence/absence data.  It is defined
-  as \eqn{1 - prob(j)}, or based on the probability of observing at
-  least \eqn{j} species in shared in compared communities.  Legendre &
-  Legendre (1998) suggest using simulations to assess the probability,
-  but the current function uses analytic result from hypergeometric
-  distribution (\code{\link{phyper}}) instead.  This probability (and
-  the index) is dependent on the number of species missing in both
+  Raup--Crick dissimilarity (\code{method = "raup"}) is a probabilistic
+  index based on presence/absence data.  It is defined as \eqn{1 -
+  prob(j)}, or based on the probability of observing at least \eqn{j}
+  species in shared in compared communities.  The current function uses
+  analytic result from hypergeometric distribution
+  (\code{\link{phyper}}) to find the probablities.  This probability
+  (and the index) is dependent on the number of species missing in both
   sites, and adding all-zero species to the data or removing missing
-  species from the data will influence the index.  The probability
-  (and the index) may be almost zero or almost one for a wide range of
+  species from the data will influence the index.  The probability (and
+  the index) may be almost zero or almost one for a wide range of
   parameter values.  The index is nonmetric: two communities with no
   shared species may have a dissimilarity slightly below one, and two
-  identical communities may have dissimilarity slightly above
-  zero. The index uses equal occurrence probabilities for all species,
-  but Raup and Crick originally suggested that sampling probabilities
-  should be proportional to species frequencies (Chase et al. 2011). A
-  simulation approach with unequal species sampling probabilities is
-  implemented in \code{\link{raupcrick}} function following Chase et
-  al. (2011).
+  identical communities may have dissimilarity slightly above zero. The
+  index uses equal occurrence probabilities for all species, but Raup
+  and Crick originally suggested that sampling probabilities should be
+  proportional to species frequencies (Chase et al. 2011). A simulation
+  approach with unequal species sampling probabilities is implemented in
+  \code{\link{raupcrick}} function following Chase et al. (2011).
   
   Chao index tries to take into account the number of unseen species
   pairs, similarly as in \code{method = "chao"} in
@@ -274,9 +272,6 @@
 
   Krebs, C. J. (1999). \emph{Ecological Methodology.} Addison Wesley Longman.
 
-  Legendre, P, & Legendre, L. (1998) \emph{Numerical Ecology}. 2nd English
-  Edition. Elsevier.
-  
   Mountford, M. D. (1962). An index of similarity and its application to
   classification problems. In: P.W.Murphy (ed.),
   \emph{Progress in Soil Zoology}, 43--50. Butterworths.

-- 
Community Ecology Package for R



More information about the debian-med-commit mailing list