[Vegan-commits] r1360 - in pkg/vegan: R inst

[noreply at r-forge.r-project.org]

Author: jarioksa
Date: 2010-11-10 11:18:55 +0100 (Wed, 10 Nov 2010)
New Revision: 1360

Modified:
   pkg/vegan/R/SSlomolino.R
   pkg/vegan/inst/ChangeLog
Log:
improved starting values for Lomolino model: now often works

Modified: pkg/vegan/R/SSlomolino.R
===================================================================
--- pkg/vegan/R/SSlomolino.R	2010-11-10 09:02:23 UTC (rev 1359)
+++ pkg/vegan/R/SSlomolino.R	2010-11-10 10:18:55 UTC (rev 1360)
@@ -3,20 +3,21 @@
               function(mCall, data, LHS)
 {
     xy <- sortedXyData(mCall[["area"]], LHS, data)
-    ## assume Asym is const*max(S)
-    .S <- max(xy[["y"]])*1.5
-    ## approximate using y = log(Smax/S - 1)
+    ## approximate with Arrhenius model on log-log
+    .p <- coef(lm(log(xy[["y"]]) ~ log(xy[["x"]])))
+    ## Asym is value at max(x) but > max(y) and xmid is x which gives
+    ## Asym/2
+    .Smax <- max(xy[["y"]])*1.1
+    .S <- exp(.p[1] + log(max(xy[["x"]])) * (.p[2]))
+    .S <- max(.S, .Smax)
+    .xmid <- exp((log(.S/2) - .p[1])/.p[2])
+    ## approximate slope using y = log(Smax/S - 1)
     .y <- log(.S/xy[["y"]] - 1)
     .x <- xy[["x"]]
-    ## xmid is the inflection point: fit a parabola to log and find it
-    ## there
-    .p <- coef(lm(.y ~ .x + I(.x^2)))
-    .xmid <- -(.p[2])/2/.p[3] - sqrt(abs(1/2/.p[3]))
-     ## estimate slope assuming Asym and xmid are known
     .z <- log(.xmid/xy[["x"]])
     .b <- coef(lm(.y ~ .z))
     ## Adjust Asym: half of y = Asym/2 at xmid
-    .S <- .S * exp(-0.5 * (.b[1]))
+    ##.S <- .S * exp(-0.5 * (.b[1]))
     value <- c(.S, .xmid, exp(.b[2]))
     names(value) <- mCall[c("Asym","xmid", "slope")]
     value

Modified: pkg/vegan/inst/ChangeLog
===================================================================
--- pkg/vegan/inst/ChangeLog	2010-11-10 09:02:23 UTC (rev 1359)
+++ pkg/vegan/inst/ChangeLog	2010-11-10 10:18:55 UTC (rev 1360)
@@ -4,6 +4,13 @@
 
 Version 1.18-16 (opened November 9, 20109
 
+	* SSlomolino: improved starting values for 'xmid' (and
+	'Asym'). Now fitspecaccum(..., "lomolino") works in several cases,
+	including BCI and bryceveg (but fails in <1% of cases). Now 'Asym'
+	and 'xmid' are estimated from Arrhenius at log-log scale: 'Asym'
+	is the predicted value at max(x), and 'xmid' is the value of x
+	giving 'Asym/2'.
+
 	* capscale: defines total inertia as the sum of absolute values of
 	eigenvalues to be consistent with cmdscale(), wcmdscale(),
 	eigenvals(), Gower and Mardia, Kent & Bibby.