[Vegan-commits] r1353 - pkg/vegan/R

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

Author: jarioksa
Date: 2010-11-07 08:25:19 +0100 (Sun, 07 Nov 2010)
New Revision: 1353

Modified:
   pkg/vegan/R/SSlomolino.R
Log:
improved starting values for SSlomollino -- still fail with fitspecaccum(BCI, "lomolino", algorith = "port")

Modified: pkg/vegan/R/SSlomolino.R
===================================================================
--- pkg/vegan/R/SSlomolino.R	2010-11-07 07:20:34 UTC (rev 1352)
+++ pkg/vegan/R/SSlomolino.R	2010-11-07 07:25:19 UTC (rev 1353)
@@ -3,12 +3,19 @@
               function(mCall, data, LHS)
 {
     xy <- sortedXyData(mCall[["area"]], LHS, data)
-    ## wild guess for starting values: assume Asym is 1.2*max(S), fit
-    ## log(S) ~ log(area), estimate xmid as area of that model giving
-    ## S/2 and assume slope is 1/slope of that model.
-    .S <- max(xy[,"y"])
-    .b <- as.vector(coef(lm(log(xy[,"y"]) ~ log(xy[,"x"]))))
-    value <- c(.S*1.2, exp(log(.S/2)/.b[2] - .b[1]), 1/.b[2])
+    ## assume Asym is const*max(S)
+    .S <- max(xy[["y"]])*1.5
+    ## approximate 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 <- exp(coef(lm(.y ~ .z - 1)))
+    value <- c(.S, .xmid, .b)
     names(value) <- mCall[c("Asym","xmid", "slope")]
     value
 },