Changed the way confidence intervals are being calculated. Now

the function confint.nls from the MASS package is being used,
after stumbling across the relevant paragraph of the MASS book.



git-svn-id: http://kriemhild.uft.uni-bremen.de/svn/drfit@76 d1b72e20-2ee0-0310-a1c4-ad5adbbefcdc

8 files changed, 181 insertions, 135 deletions

diff --git a/DESCRIPTION b/DESCRIPTION
index 5d92da7..e12b6bf 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,14 +1,14 @@
 Package: drfit
-Version: 0.05-75
-Date: 2006-04-28
+Version: 0.05-76
+Date: 2006-05-04
 Title: Dose-response data evaluation
 Author: Johannes Ranke <jranke@uni-bremen.de>
 Maintainer: Johannes Ranke <jranke@uni-bremen.de>
-Depends: R (>= 2.1.0),stats,RODBC
+Depends: R (>= 2.1.0),stats,MASS,RODBC
 Description: drfit provides basic and easy-to-use functions for fitting
 	dose-response curves to dose-response data, calculating some
 	(eco)toxicological parameters and plotting the results. Functions that are
-	fitted are the cumulative densitiy function of the lognormal distribution
+	fitted are the cumulative density function of the lognormal distribution
 	(probit fit), of the logistic distribution (logit fit), of the weibull
 	distribution (weibull fit) and a linear-logistic model ("linlogit" fit),
 	derived from the latter, which is used to describe data showing stimulation
diff --git a/R/drfit.R b/R/drfit.R
index 1bbffa8..b33d754 100644
--- a/R/drfit.R
+++ b/R/drfit.R
@@ -1,27 +1,28 @@
 drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
         probit = TRUE, logit = FALSE, weibull = FALSE,
-        linlogit = FALSE, conf = FALSE,
+        linlogit = FALSE, level = 0.95,
         linlogitWrong = NA, allWrong = NA,
         s0 = 0.5, b0 = 2, f0 = 0)
 {
-    if(!is.null(data$ok)) data <- subset(data,ok!="no fit") # Don't use data with 
-                                                            # ok set to "no fit"
+    require(MASS)
+    if(!is.null(data$ok)) data <- subset(data,ok!="no fit") # Don't use data
+                                                            # with ok set to
+                                                            # "no fit"
     substances <- levels(data$substance)
 
-    ri <- rix <- 0                        # ri is the index over the result rows
-                                          # rix is used later to check if any
-                                          # model result was appended
-    rsubstance <- array()                 # the substance names in the results
-    rndl <- vector()                      # number of dose levels
-    rn <- vector()                        # mean number of replicates 
-                                          # in each dose level
-    runit <- vector()                     # vector of units for each result row
-    rlhd <- rlld <- vector()              # highest and lowest doses tested
-    mtype <- array()                      # the modeltypes
-    sigma <- array()                      # the standard deviation of the residuals
+    ri <- rix <- 0                  # ri is the index over the result rows
+                                    # rix is used later to check if any
+                                    # model result was appended
+    rsubstance <- array()           # the substance names in the results
+    rndl <- vector()                # number of dose levels
+    rn <- vector()                  # mean number of replicates 
+                                    # in each dose level
+    runit <- vector()               # vector of units for each result row
+    rlhd <- rlld <- vector()        # highest and lowest doses tested
+    mtype <- array()                # the modeltypes
+    sigma <- array()                # the standard deviation of the residuals
     logED50 <- vector()
-    stderrlogED50 <- vector()
-    conflogED50 <- vector()
+    logED50low <- logED50high <- vector()
     a <- b <- c <- vector()
 
     splitted <- split(data,data$substance)
@@ -87,15 +88,21 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
                             if (logED50[[ri]] > rlhd[[ri]]) {
                                 mtype[[ri]] <- "no fit"
                                 logED50[[ri]] <- NA
-                                stderrlogED50[[ri]] <- NA
-                                conflogED50[[ri]] <- NA
+                                logED50low[[ri]] <- NA
+                                logED50high[[ri]] <- NA
                                 a[[ri]] <- NA
                                 b[[ri]] <- NA
                                 c[[ri]] <- NA
                             } else {
                                 mtype[[ri]] <- "linlogit"
-                                stderrlogED50[[ri]] <- s$parameters["logED50","Std. Error"]
-                                conflogED50[[ri]] <- stderrlogED50[[ri]] * qt(0.975, n - 3)
+                                logED50conf <- try(confint(m,"logED50",level=level))
+                                if (!inherits(logED50conf, "try-error")) {
+                                    logED50low[[ri]] <- logED50conf[[1]]
+                                    logED50high[[ri]] <- logED50conf[[2]]
+                                } else {
+                                    logED50low[[ri]] <- NA
+                                    logED50high[[ri]] <- NA
+                                }
                                 a[[ri]] <- coef(m)[["logED50"]]
                                 b[[ri]] <- coef(m)[["b"]]
                                 c[[ri]] <- coef(m)[["f"]]
@@ -125,14 +132,20 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
                                 if (logED50[[ri]] > rlhd[[ri]]) {
                                     mtype[[ri]] <- "no fit"
                                     logED50[[ri]] <- NA
-                                    stderrlogED50[[ri]] <- NA
-                                    conflogED50[[ri]] <- NA
+                                    logED50low[[ri]] <- NA
+                                    logED50high[[ri]] <- NA
                                     a[[ri]] <- NA
                                     b[[ri]] <- NA
                                 } else {
                                     mtype[[ri]] <- "probit"
-                                    stderrlogED50[[ri]] <- s$parameters["logED50","Std. Error"]
-                                    conflogED50[[ri]] <- stderrlogED50[[ri]] * qt(0.975, n - 2)
+                                    logED50conf <- try(confint(m,"logED50",level=level))
+                                    if (!inherits(logED50conf, "try-error")) {
+                                        logED50low[[ri]] <- logED50conf[[1]]
+                                        logED50high[[ri]] <- logED50conf[[2]]
+                                    } else {
+                                        logED50low[[ri]] <- NA
+                                        logED50high[[ri]] <- NA
+                                    }
                                     a[[ri]] <- coef(m)[["logED50"]]
                                     b[[ri]] <- coef(m)[["scale"]]
                                 }
@@ -163,14 +176,20 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
                                 if (logED50[[ri]] > rlhd[[ri]]) {
                                     mtype[[ri]] <- "no fit"
                                     logED50[[ri]] <- NA
-                                    stderrlogED50[[ri]] <- NA
-                                    conflogED50[[ri]] <- NA
+                                    logED50low[[ri]] <- NA
+                                    logED50high[[ri]] <- NA
                                     a[[ri]] <- NA
                                     b[[ri]] <- NA
                                 } else {
                                     mtype[[ri]] <- "logit"
-                                    stderrlogED50[[ri]] <- s$parameters["logED50","Std. Error"]
-                                    conflogED50[[ri]] <- stderrlogED50[[ri]] * qt(0.975, n - 2)
+                                    logED50conf <- try(confint(m,"logED50",level=level))
+                                    if (!inherits(logED50conf, "try-error")) {
+                                        logED50low[[ri]] <- logED50conf[[1]]
+                                        logED50high[[ri]] <- logED50conf[[2]]
+                                    } else {
+                                        logED50low[[ri]] <- NA
+                                        logED50high[[ri]] <- NA
+                                    }
                                 }
                             }
                         }
@@ -200,17 +219,15 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
                                 }
                                 logED50[[ri]] <- nlm(sqrdev,startlogED50[[i]])$estimate
                                 c[[ri]] <- NA
+                                logED50low[[ri]] <- NA
+                                logED50high[[ri]] <- NA
                                 if (logED50[[ri]] > rlhd[[ri]]) {
                                     mtype[[ri]] <- "no fit"
                                     logED50[[ri]] <- NA
-                                    stderrlogED50[[ri]] <- NA
-                                    conflogED50[[ri]] <- NA
                                     a[[ri]] <- NA
                                     b[[ri]] <- NA
                                 } else {
                                     mtype[[ri]] <- "weibull"
-                                    stderrlogED50[[ri]] <- NA
-                                    conflogED50[[ri]] <- stderrlogED50[[ri]] * qt(0.975, n - 2)
                                 }
                             }
                         }
@@ -247,21 +264,17 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
             }
             sigma[[ri]] <- NA
             logED50[[ri]] <- NA
-            stderrlogED50[[ri]] <- NA
-            conflogED50[[ri]] <- NA
+            logED50low[[ri]] <- NA
+            logED50high[[ri]] <- NA
             a[[ri]] <- NA
             b[[ri]] <- NA
             c[[ri]] <- NA
         }
     }
-    if (conf)
-    {
-        results <- data.frame(rsubstance, rndl, rn, rlld, rlhd, mtype, logED50, conflogED50, runit, sigma, a, b)
-        names(results) <- c("Substance","ndl","n","lld","lhd","mtype","logED50","conf","unit","sigma","a","b")
-    } else {
-        results <- data.frame(rsubstance, rndl, rn, rlld, rlhd, mtype, logED50, stderrlogED50, runit, sigma, a, b)
-        names(results) <- c("Substance","ndl","n","lld","lhd","mtype","logED50","std","unit","sigma","a","b")
-    } 
+    results <- data.frame(rsubstance, rndl, rn, rlld, rlhd, mtype, 
+        logED50, logED50low, logED50high, runit, sigma, a, b)
+    names(results) <- c("Substance","ndl","n","lld","lhd","mtype",
+        "logED50",names(logED50conf)[[1]],names(logED50conf)[[2]],"unit","sigma","a","b")
 
     if (linlogit) {
         results$c <- c
diff --git a/man/drfit.Rd b/man/drfit.Rd
index a89b5dd..453c0f8 100644
--- a/man/drfit.Rd
+++ b/man/drfit.Rd
@@ -7,7 +7,7 @@
 }
 \usage{
   drfit(data, startlogED50 = NA, chooseone = TRUE, probit = TRUE, logit = FALSE,
-  weibull = FALSE, linlogit = FALSE, conf = FALSE, linlogitWrong = NA,
+  weibull = FALSE, linlogit = FALSE, level = 0.95, linlogitWrong = NA,
   allWrong = NA, s0 = 0.5, b0 = 2, f0 = 0)
 }
 \arguments{
@@ -43,10 +43,8 @@
     A boolean defining if the linear-logistic function
     \code{\link{linlogitf}} as defined by van Ewijk and Hoekstra 1993 is
     fitted to the data. Default is FALSE.}
-  \item{conf}{
-    A boolean defining if a confidence interval of the log ED50 is to be
-    listed for alpha = 5 \% (two-sided). If FALSE (default), the standard
-    deviation is listed in the output of drfit.}
+  \item{level}{
+    The level for the confidence interval listed for the log ED50.}
   \item{linlogitWrong}{
     An optional vector containing the names of the substances for which the
     linlogit function produces a wrong fit.}
@@ -81,10 +79,10 @@
     \code{lhd} is the decadic logarithm of the highest dose.  For the
     \dQuote{linlogit}, \dQuote{logit} and \dQuote{probit} models, the
     parameter \code{a} that is reported coincides with the logED50, i.e the
-    logED50 is one of the model parameters that is being fitted, and
-    therefore, depending on the argument \code{conf}, a standard deviation
-    \code{std} or a confidence interval \code{conf} is reported for the
-    logED50. In the case of the \dQuote{weibull} model, \code{a} is a
+    logED50 is one of the model parameters that is being fitted. Therefore,
+    a confidence interval for the confidence level \code{level} is calculated
+    using the \code{\link[MASS:confint]{confint.nls}} function and listed.
+    In the case of the \dQuote{weibull} model, \code{a} is a
     location parameter.  Parameter \code{b} in the case of the
     \dQuote{linlogit} fit is the variable b from the \code{\link{linlogitf}}
     function. In the case of \dQuote{probit} fit it is the standard deviation
diff --git a/tests/IM1xIPC81.Rout.save b/tests/IM1xIPC81.Rout.save
index b67ed2a..b56b8c9 100644
--- a/tests/IM1xIPC81.Rout.save
+++ b/tests/IM1xIPC81.Rout.save
@@ -16,6 +16,7 @@ Type 'demo()' for some demos, 'help()' for on-line help, or
 Type 'q()' to quit R.
 
 > library(drfit)
+Loading required package: MASS
 Loading required package: RODBC
 > data(IM1xIPC81)
 > rIM1xIPC81 <- drfit(IM1xIPC81,linlogit=TRUE)
@@ -35,41 +36,45 @@ IM15 BF4: Fitting data...
 IM16 BF4: Fitting data...
 
 IM17 BF4: Fitting data...
+Waiting for profiling to be done...
 Error in numericDeriv(form[[3]], names(ind), env) : 
 	Missing value or an infinity produced when evaluating the model
 In addition: Warning message:
 NaNs produced in: pnorm(q, mean, sd, lower.tail, log.p) 
 
 IM18 BF4: Fitting data...
+Waiting for profiling to be done...
 Error in numericDeriv(form[[3]], names(ind), env) : 
 	Missing value or an infinity produced when evaluating the model
 In addition: Warning message:
 NaNs produced in: pnorm(q, mean, sd, lower.tail, log.p) 
 
 IM19 BF4: Fitting data...
+Waiting for profiling to be done...
 Error in numericDeriv(form[[3]], names(ind), env) : 
 	Missing value or an infinity produced when evaluating the model
 In addition: Warning message:
 NaNs produced in: pnorm(q, mean, sd, lower.tail, log.p) 
 
 IM1-10 BF4: Fitting data...
+Waiting for profiling to be done...
 > rIM1xIPC81
-   Substance ndl   n         lld      lhd    mtype   logED50        std unit
-1   IM13 BF4   9  81  0.59176003 3.000000 inactive        NA         NA   然
-2   IM14 BF4  20 216 -0.01030017 3.176091   no fit        NA         NA   然
-3   IM15 BF4   9 135  0.59176003 3.000000 inactive        NA         NA   然
-4   IM16 BF4   9 108  0.59176003 3.000000 inactive        NA         NA   然
-5   IM17 BF4   9  81  0.59176003 3.000000 linlogit 2.5785719 0.03819482   然
-6   IM18 BF4   9 135  0.59176003 3.000000 linlogit 1.6806040 0.02982682   然
-7   IM19 BF4   9  81  0.59176003 3.000000 linlogit 1.6496121 0.02636274   然
-8 IM1-10 BF4  11 162 -0.01030017 3.000000 linlogit 0.7696961 0.04165397   然
-      sigma         a        b          c
-1        NA        NA       NA         NA
-2        NA        NA       NA         NA
-3        NA        NA       NA         NA
-4        NA        NA       NA         NA
-5 0.2375857 2.5785719 2.300216 0.01468034
-6 0.2324871 1.6806040 2.237338 0.05719295
-7 0.1452716 1.6496121 1.976976 0.10956088
-8 0.2987708 0.7696961 1.936253 0.45808611
+   Substance ndl   n         lld      lhd    mtype   logED50      2.5%
+1   IM13 BF4   9  81  0.59176003 3.000000 inactive        NA        NA
+2   IM14 BF4  20 216 -0.01030017 3.176091   no fit        NA        NA
+3   IM15 BF4   9 135  0.59176003 3.000000 inactive        NA        NA
+4   IM16 BF4   9 108  0.59176003 3.000000 inactive        NA        NA
+5   IM17 BF4   9  81  0.59176003 3.000000 linlogit 2.5785719 2.5057219
+6   IM18 BF4   9 135  0.59176003 3.000000 linlogit 1.6806040 1.6229756
+7   IM19 BF4   9  81  0.59176003 3.000000 linlogit 1.6496121 1.5984793
+8 IM1-10 BF4  11 162 -0.01030017 3.000000 linlogit 0.7696961 0.6870081
+      97.5% unit     sigma         a        b          c
+1        NA   然        NA        NA       NA         NA
+2        NA   然        NA        NA       NA         NA
+3        NA   然        NA        NA       NA         NA
+4        NA   然        NA        NA       NA         NA
+5 2.6617385   然 0.2375857 2.5785719 2.300216 0.01468034
+6 1.7418917   然 0.2324871 1.6806040 2.237338 0.05719295
+7 1.7031463   然 0.1452716 1.6496121 1.976976 0.10956088
+8 0.8544213   然 0.2987708 0.7696961 1.936253 0.45808612
 > 
diff --git a/tests/IM1xVibrio.Rout.save b/tests/IM1xVibrio.Rout.save
index 5d37703..bebdcea 100644
--- a/tests/IM1xVibrio.Rout.save
+++ b/tests/IM1xVibrio.Rout.save
@@ -16,42 +16,51 @@ Type 'demo()' for some demos, 'help()' for on-line help, or
 Type 'q()' to quit R.
 
 > library(drfit)
+Loading required package: MASS
 Loading required package: RODBC
 > data(IM1xVibrio)
 > rIM1xVibrio <- drfit(IM1xVibrio)
 
 IM13 BF4: Fitting data...
+Waiting for profiling to be done...
 
 IM14 BF4: Fitting data...
+Waiting for profiling to be done...
 
 IM15 BF4: Fitting data...
+Waiting for profiling to be done...
 
 IM16 BF4: Fitting data...
+Waiting for profiling to be done...
 
 IM17 BF4: Fitting data...
+Waiting for profiling to be done...
 
 IM18 BF4: Fitting data...
+Waiting for profiling to be done...
 
 IM19 BF4: Fitting data...
+Waiting for profiling to be done...
 
 IM1-10 BF4: Fitting data...
+Waiting for profiling to be done...
 > rIM1xVibrio
-   Substance ndl  n      lld     lhd  mtype    logED50        std unit
-1   IM13 BF4  13 22 -5.30103 4.69897 probit  3.9399128 0.03658457   然
-2   IM14 BF4  12 57 -4.30103 4.69897 probit  3.5442110 0.02060194   然
-3   IM15 BF4   7 20 -0.30103 4.69897 probit  3.1398488 0.01254468   然
-4   IM16 BF4  13 30 -4.30103 4.69897 probit  3.1804247 0.01801591   然
-5   IM17 BF4  10 20 -1.30103 4.69897 probit  2.4317301 0.03015341   然
-6   IM18 BF4  10 17 -2.30103 2.69897 probit  1.4015187 0.04053083   然
-7   IM19 BF4  13 22 -5.30103 2.69897 probit  0.7157942 0.02085130   然
-8 IM1-10 BF4  11 20 -6.30103 2.69897 probit -0.1790307 0.03623821   然
-       sigma          a         b
-1 0.06302893  3.9399128 0.5763074
-2 0.05905375  3.5442110 0.6561300
-3 0.02695476  3.1398488 0.5526675
-4 0.04514079  3.1804247 0.4881559
-5 0.04773317  2.4317301 0.6386547
-6 0.05611347  1.4015187 0.8008513
-7 0.02956284  0.7157942 0.6969275
-8 0.04098662 -0.1790307 0.8019351
+   Substance ndl  n      lld     lhd  mtype    logED50       2.5%      97.5%
+1   IM13 BF4  13 22 -5.30103 4.69897 probit  3.9399128  3.8636215  4.0160099
+2   IM14 BF4  12 57 -4.30103 4.69897 probit  3.5442110  3.5030166  3.5855474
+3   IM15 BF4   7 20 -0.30103 4.69897 probit  3.1398488  3.1135935  3.1662394
+4   IM16 BF4  13 30 -4.30103 4.69897 probit  3.1804247  3.1435815  3.2176152
+5   IM17 BF4  10 20 -1.30103 4.69897 probit  2.4317301  2.3680656  2.4951179
+6   IM18 BF4  10 17 -2.30103 2.69897 probit  1.4015187  1.3146736  1.4880647
+7   IM19 BF4  13 22 -5.30103 2.69897 probit  0.7157942  0.6703743  0.7592458
+8 IM1-10 BF4  11 20 -6.30103 2.69897 probit -0.1790307 -0.2569430 -0.1037183
+  unit      sigma          a         b
+1   然 0.06302893  3.9399128 0.5763074
+2   然 0.05905375  3.5442110 0.6561300
+3   然 0.02695476  3.1398488 0.5526675
+4   然 0.04514079  3.1804247 0.4881559
+5   然 0.04773317  2.4317301 0.6386547
+6   然 0.05611347  1.4015187 0.8008513
+7   然 0.02956284  0.7157942 0.6969275
+8   然 0.04098662 -0.1790307 0.8019351
 > 
diff --git a/tests/XY.Rout.save b/tests/XY.Rout.save
index 7bcc0df..d8098ae 100644
--- a/tests/XY.Rout.save
+++ b/tests/XY.Rout.save
@@ -16,6 +16,7 @@ Type 'demo()' for some demos, 'help()' for on-line help, or
 Type 'q()' to quit R.
 
 > library(drfit)
+Loading required package: MASS
 Loading required package: RODBC
 > data(XY)
 > rXY <- drfit(XY)
@@ -23,15 +24,16 @@ Loading required package: RODBC
 Control: Fitting data...
 
 Substance X: Fitting data...
+Waiting for profiling to be done...
 
 Substance Y: Fitting data...
 > rXY
-    Substance ndl  n  lld  lhd    mtype  logED50        std unit      sigma
-1     Control   1  6 -Inf -Inf inactive       NA         NA mg/L         NA
-2 Substance X   4 12    1    3   probit 2.161905 0.02425144 mg/L 0.04131746
-3 Substance Y   4 12    1    3   active       NA         NA mg/L         NA
-         a         b
-1       NA        NA
-2 2.161905 0.5061306
-3       NA        NA
+    Substance ndl  n  lld  lhd    mtype  logED50     2.5%    97.5% unit
+1     Control   1  6 -Inf -Inf inactive       NA       NA       NA mg/L
+2 Substance X   4 12    1    3   probit 2.161905 2.105861 2.214359 mg/L
+3 Substance Y   4 12    1    3   active       NA       NA       NA mg/L
+       sigma        a         b
+1         NA       NA        NA
+2 0.04131746 2.161905 0.5061306
+3         NA       NA        NA
 > 
diff --git a/tests/antifoul.Rout.save b/tests/antifoul.Rout.save
index 5bc166c..9967240 100644
--- a/tests/antifoul.Rout.save
+++ b/tests/antifoul.Rout.save
@@ -16,18 +16,21 @@ Type 'demo()' for some demos, 'help()' for on-line help, or
 Type 'q()' to quit R.
 
 > library(drfit)
+Loading required package: MASS
 Loading required package: RODBC
 > data(antifoul)
 > rantifoul <- drfit(antifoul)
 
 TBT: Fitting data...
+Waiting for profiling to be done...
 
 Zn Pyrithion: Fitting data...
+Waiting for profiling to be done...
 > rantifoul
-     Substance ndl   n       lld     lhd  mtype    logED50        std unit
-1          TBT  38 135 -2.709271 2.39794 probit -0.1643560 0.05349216   然
-2 Zn Pyrithion  27  81 -2.107210 2.00000 probit -0.3979222 0.05393693   然
-      sigma          a         b
-1 0.1928611 -0.1643560 0.6763882
-2 0.2286560 -0.3979222 0.4169632
+     Substance ndl   n       lld     lhd  mtype    logED50       2.5%
+1          TBT  38 135 -2.709271 2.39794 probit -0.1643560 -0.2650411
+2 Zn Pyrithion  27  81 -2.107210 2.00000 probit -0.3979222 -0.5078291
+        97.5% unit     sigma          a         b
+1 -0.05648615   然 0.1928611 -0.1643560 0.6763882
+2 -0.29156056   然 0.2286560 -0.3979222 0.4169632
 > 
diff --git a/tests/pyrithione.Rout.save b/tests/pyrithione.Rout.save
index fb11ca9..2a05dd8 100644
--- a/tests/pyrithione.Rout.save
+++ b/tests/pyrithione.Rout.save
@@ -16,70 +16,86 @@ Type 'demo()' for some demos, 'help()' for on-line help, or
 Type 'q()' to quit R.
 
 > library(drfit)
+Loading required package: MASS
 Loading required package: RODBC
 > data(pyrithione)
 > rpyr <- drfit(pyrithione,linlogit=TRUE,linlogitWrong=c("MSPT","MSPHI"))
 
 Na Pyrithion: Fitting data...
+Waiting for profiling to be done...
+Error in prof$getProfile() : number of iterations exceeded maximum of 50
 
 Pyridin: Fitting data...
 
 PyNO: Fitting data...
 
 (PT)2: Fitting data...
+Waiting for profiling to be done...
 
 MSO2P: Fitting data...
 
 MSPHI: Fitting data...
+Waiting for profiling to be done...
 
 PyS: Fitting data...
 Error in nls(response ~ linlogitf(dose, 1, f, logED50, b), data = tmp,  : 
 	step factor 0.000488281 reduced below 'minFactor' of 0.000976562
+Waiting for profiling to be done...
+Error in numericDeriv(form[[3]], names(ind), env) : 
+	Missing value or an infinity produced when evaluating the model
+In addition: Warning message:
+NaNs produced in: pnorm(q, mean, sd, lower.tail, log.p) 
 
 Zn Pyrithion: Fitting data...
+Waiting for profiling to be done...
 
 Cu Pyrithion: Fitting data...
+Waiting for profiling to be done...
 
 Fe Pyrithion: Fitting data...
+Waiting for profiling to be done...
 
 (PyS)2: Fitting data...
 Error in nls(response ~ linlogitf(dose, 1, f, logED50, b), data = tmp,  : 
 	singular gradient
+Waiting for profiling to be done...
 
 MSPT: Fitting data...
+Waiting for profiling to be done...
 
 TBT: Fitting data...
+Waiting for profiling to be done...
 
 NaJ: Fitting data...
 > rpyr
-      Substance ndl   n       lld     lhd    mtype    logED50        std unit
-1  Na Pyrithion  20 108 -2.107210 2.00000 linlogit -0.3461318 0.04237833   然
-2       Pyridin  19 161  0.591760 3.69897 inactive         NA         NA   然
-3          PyNO  19  81  0.591760 3.69897   no fit         NA         NA   然
-4         (PT)2  19  81 -2.408240 2.00000 linlogit -0.4197031 0.05030460   然
-5         MSO2P  28 108 -1.709271 3.69897 inactive         NA         NA   然
-6         MSPHI  19  81  0.591760 3.69897   probit  3.6625025 0.06724673   然
-7           PyS  18  80 -0.408240 3.00000   probit  2.6817944 0.02031252   然
-8  Zn Pyrithion  27  81 -2.107210 2.00000 linlogit -0.4132885 0.06234927   然
-9  Cu Pyrithion  19  79 -2.408240 2.00000 linlogit -0.3074717 0.07918130   然
-10 Fe Pyrithion  19  81 -2.408240 2.00000 linlogit -0.3534572 0.07263386   然
-11       (PyS)2  20  81 -0.408240 3.00000   probit  1.7840499 0.05230663   然
-12         MSPT  18 108 -0.408240 3.00000   probit  2.1384811 0.03379904   然
-13          TBT  38 135 -2.709271 2.39794 linlogit -0.1580954 0.05474371   然
-14          NaJ  10 108  0.591760 3.30103 inactive         NA         NA   然
-        sigma          a         b           c
-1  0.20857751 -0.3461318 2.1188272 -0.26698979
-2          NA         NA        NA          NA
-3  0.09188853         NA        NA          NA
-4  0.20509070 -0.4197031 1.8466994 -0.51508088
-5          NA         NA        NA          NA
-6  0.12943273  3.6625025 0.4387540          NA
-7  0.13401920  2.6817944 0.1250523          NA
-8  0.22877916 -0.4132885 1.7424846  0.61078660
-9  0.24494451 -0.3074717 1.0150649 -0.04897422
-10 0.23383185 -0.3534572 1.1823558  0.02510613
-11 0.23298721  1.7840499 0.4098018          NA
-12 0.22881257  2.1384811 0.2112625          NA
-13 0.19264079 -0.1580954 1.0406576 -0.02416104
-14         NA         NA        NA          NA
+      Substance ndl   n       lld     lhd    mtype    logED50       2.5%
+1  Na Pyrithion  20 108 -2.107210 2.00000 linlogit -0.3461318         NA
+2       Pyridin  19 161  0.591760 3.69897 inactive         NA         NA
+3          PyNO  19  81  0.591760 3.69897   no fit         NA         NA
+4         (PT)2  19  81 -2.408240 2.00000 linlogit -0.4197042 -0.5326102
+5         MSO2P  28 108 -1.709271 3.69897 inactive         NA         NA
+6         MSPHI  19  81  0.591760 3.69897   probit  3.6625025  3.5548278
+7           PyS  18  80 -0.408240 3.00000   probit  2.6817944         NA
+8  Zn Pyrithion  27  81 -2.107210 2.00000 linlogit -0.4132886 -0.5293327
+9  Cu Pyrithion  19  79 -2.408240 2.00000 linlogit -0.3074717 -0.4753489
+10 Fe Pyrithion  19  81 -2.408240 2.00000 linlogit -0.3534573 -0.5003829
+11       (PyS)2  20  81 -0.408240 3.00000   probit  1.7840499  1.6726792
+12         MSPT  18 108 -0.408240 3.00000   probit  2.1384811  2.0709306
+13          TBT  38 135 -2.709271 2.39794 linlogit -0.1580954 -0.2602139
+14          NaJ  10 108  0.591760 3.30103 inactive         NA         NA
+         97.5% unit      sigma          a         b           c
+1           NA   然 0.20857751 -0.3461318 2.1188275 -0.26698926
+2           NA   然         NA         NA        NA          NA
+3           NA   然 0.09188853         NA        NA          NA
+4           NA   然 0.20509070 -0.4197042 1.8467290 -0.51505615
+5           NA   然         NA         NA        NA          NA
+6   3.88053203   然 0.12943273  3.6625025 0.4387540          NA
+7           NA   然 0.13401920  2.6817944 0.1250523          NA
+8  -0.29979846   然 0.22877916 -0.4132886 1.7424848  0.61079006
+9  -0.15442805   然 0.24494451 -0.3074717 1.0150650 -0.04897413
+10 -0.21459555   然 0.23383185 -0.3534573 1.1823561  0.02510663
+11  1.88824519   然 0.23298721  1.7840499 0.4098018          NA
+12  2.20511512   然 0.22881257  2.1384811 0.2112625          NA
+13 -0.05045493   然 0.19264079 -0.1580954 1.0406577 -0.02416082
+14          NA   然         NA         NA        NA          NA
 >