Changed the terminology from EC50 to ED50 reflecting that the package is meant

for dose-response analysis in general and not only for concentration-response
analysis. This is the first revision that has been prepared with the knowledge
of the existence of the drc package for dose-response curve analysis.



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

1 files changed, 44 insertions, 44 deletions

diff --git a/R/drfit.R b/R/drfit.R
index d8631e6..e74667d 100644
--- a/R/drfit.R
+++ b/R/drfit.R
@@ -44,7 +44,7 @@ linlogitf <- function(x,k,f,mu,b)
     k*(1 + f*x) / (1 + ((2*f*(10^mu) + 1) * ((x/(10^mu))^b)))
 }
 
-drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
+drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
         probit = TRUE, logit = FALSE, weibull = FALSE,
         linlogit = FALSE, linlogitWrong = NA, allWrong = NA,
         s0 = 0.5, b0 = 2, f0 = 0)
@@ -62,8 +62,8 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
     rlhd <- rlld <- vector()              # highest and lowest doses tested
     mtype <- array()                      # the modeltypes
     sigma <- array()                      # the standard deviation of the residuals
-    logEC50 <- vector()
-    stderrlogEC50 <- vector()
+    logED50 <- vector()
+    stderrlogED50 <- vector()
     a <- b <- c <- vector()
 
     splitted <- split(data,data$substance)
@@ -89,9 +89,9 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
         }
         rix <- ri
         if (!nodata) {
-            if (is.na(startlogEC50[i])){
+            if (is.na(startlogED50[i])){
                 w <- 1/abs(tmp$response - 0.3)
-                startlogEC50[[i]] <- sum(w * log10(tmp$dose))/sum(w)
+                startlogED50[[i]] <- sum(w * log10(tmp$dose))/sum(w)
             }
             highestdose <- max(tmp$dose)
             lowestdose <- min(tmp$dose)
@@ -104,9 +104,9 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
                 if (linlogit && 
                     length(subset(linlogitWrong,linlogitWrong == i))==0 &&
                     length(subset(allWrong,allWrong == i))==0) {
-                    m <- try(nls(response ~ linlogitf(dose,1,f,logEC50,b),
+                    m <- try(nls(response ~ linlogitf(dose,1,f,logED50,b),
                             data=tmp,
-                            start=list(f=f0,logEC50=startlogEC50[[i]],b=b0)))
+                            start=list(f=f0,logED50=startlogED50[[i]],b=b0)))
                     if (!inherits(m, "try-error")) {
                         fit <- TRUE
                         ri <- ri + 1
@@ -117,18 +117,18 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
                         runit[[ri]] <- unit
                         rlld[[ri]] <- log10(lowestdose)
                         rlhd[[ri]] <- log10(highestdose)
-                        logEC50[[ri]] <- coef(m)[["logEC50"]]
-                        if (logEC50[[ri]] > rlhd[[ri]]) {
+                        logED50[[ri]] <- coef(m)[["logED50"]]
+                        if (logED50[[ri]] > rlhd[[ri]]) {
                             mtype[[ri]] <- "no fit"
-                            logEC50[[ri]] <- NA
-                            stderrlogEC50[[ri]] <- NA
+                            logED50[[ri]] <- NA
+                            stderrlogED50[[ri]] <- NA
                             a[[ri]] <- NA
                             b[[ri]] <- NA
                             c[[ri]] <- NA
                         } else {
                             mtype[[ri]] <- "linlogit"
-                            stderrlogEC50[[ri]] <- s$parameters["logEC50","Std. Error"]
-                            a[[ri]] <- coef(m)[["logEC50"]]
+                            stderrlogED50[[ri]] <- s$parameters["logED50","Std. Error"]
+                            a[[ri]] <- coef(m)[["logED50"]]
                             b[[ri]] <- coef(m)[["b"]]
                             c[[ri]] <- coef(m)[["f"]]
                         }
@@ -137,9 +137,9 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
 
                 if (probit &&
                     length(subset(allWrong,allWrong == i))==0) {
-                    m <- try(nls(response ~ pnorm(-log10(dose),-logEC50,scale),
+                    m <- try(nls(response ~ pnorm(-log10(dose),-logED50,scale),
                                 data=tmp,
-                                start=list(logEC50=startlogEC50[[i]],scale=1)))
+                                start=list(logED50=startlogED50[[i]],scale=1)))
                     if (chooseone==FALSE || fit==FALSE) {
                         if (!inherits(m, "try-error")) {
                             fit <- TRUE
@@ -151,18 +151,18 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
                             runit[[ri]] <- unit
                             rlld[[ri]] <- log10(lowestdose)
                             rlhd[[ri]] <- log10(highestdose)
-                            logEC50[[ri]] <- coef(m)[["logEC50"]]
+                            logED50[[ri]] <- coef(m)[["logED50"]]
                             c[[ri]] <- NA
-                            if (logEC50[[ri]] > rlhd[[ri]]) {
+                            if (logED50[[ri]] > rlhd[[ri]]) {
                                 mtype[[ri]] <- "no fit"
-                                logEC50[[ri]] <- NA
-                                stderrlogEC50[[ri]] <- NA
+                                logED50[[ri]] <- NA
+                                stderrlogED50[[ri]] <- NA
                                 a[[ri]] <- NA
                                 b[[ri]] <- NA
                             } else {
                                 mtype[[ri]] <- "probit"
-                                stderrlogEC50[[ri]] <- s$parameters["logEC50","Std. Error"]
-                                a[[ri]] <- coef(m)[["logEC50"]]
+                                stderrlogED50[[ri]] <- s$parameters["logED50","Std. Error"]
+                                a[[ri]] <- coef(m)[["logED50"]]
                                 b[[ri]] <- coef(m)[["scale"]]
                             }
                         }
@@ -171,9 +171,9 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
 
                 if (logit &&
                     length(subset(allWrong,allWrong == i))==0) {
-                    m <- try(nls(response ~ plogis(-log10(dose),-logEC50,scale),
+                    m <- try(nls(response ~ plogis(-log10(dose),-logED50,scale),
                             data=tmp,
-                            start=list(logEC50=startlogEC50[[i]],scale=1)))
+                            start=list(logED50=startlogED50[[i]],scale=1)))
                     if (chooseone==FALSE || fit==FALSE) {
                         if (!inherits(m, "try-error")) {
                             fit <- TRUE
@@ -185,18 +185,18 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
                             runit[[ri]] <- unit
                             rlld[[ri]] <- log10(lowestdose)
                             rlhd[[ri]] <- log10(highestdose)
-                            logEC50[[ri]] <- a[[ri]] <- coef(m)[["logEC50"]]
+                            logED50[[ri]] <- a[[ri]] <- coef(m)[["logED50"]]
                             b[[ri]] <- coef(m)[["scale"]]
                             c[[ri]] <- NA
-                            if (logEC50[[ri]] > rlhd[[ri]]) {
+                            if (logED50[[ri]] > rlhd[[ri]]) {
                                 mtype[[ri]] <- "no fit"
-                                logEC50[[ri]] <- NA
-                                stderrlogEC50[[ri]] <- NA
+                                logED50[[ri]] <- NA
+                                stderrlogED50[[ri]] <- NA
                                 a[[ri]] <- NA
                                 b[[ri]] <- NA
                             } else {
                                 mtype[[ri]] <- "logit"
-                                stderrlogEC50[[ri]] <- s$parameters["logEC50","Std. Error"]
+                                stderrlogED50[[ri]] <- s$parameters["logED50","Std. Error"]
                             }
                         }
                     }
@@ -206,7 +206,7 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
                     length(subset(allWrong,allWrong == i))==0) {
                     m <- try(nls(response ~ pweibull(-log10(dose)+location,shape),
                             data=tmp,
-                            start=list(location=startlogEC50[[i]],shape=s0)))
+                            start=list(location=startlogED50[[i]],shape=s0)))
                     if (chooseone==FALSE || fit==FALSE) {
                         if (!inherits(m, "try-error")) {
                             fit <- TRUE
@@ -223,17 +223,17 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
                             sqrdev <- function(logdose) {
                                 (0.5 - pweibull( - logdose + a[[ri]], b[[ri]]))^2 
                             }
-                            logEC50[[ri]] <- nlm(sqrdev,startlogEC50[[i]])$estimate
+                            logED50[[ri]] <- nlm(sqrdev,startlogED50[[i]])$estimate
                             c[[ri]] <- NA
-                            if (logEC50[[ri]] > rlhd[[ri]]) {
+                            if (logED50[[ri]] > rlhd[[ri]]) {
                                 mtype[[ri]] <- "no fit"
-                                logEC50[[ri]] <- NA
-                                stderrlogEC50[[ri]] <- NA
+                                logED50[[ri]] <- NA
+                                stderrlogED50[[ri]] <- NA
                                 a[[ri]] <- NA
                                 b[[ri]] <- NA
                             } else {
                                 mtype[[ri]] <- "weibull"
-                                stderrlogEC50[[ri]] <- NA
+                                stderrlogED50[[ri]] <- NA
                             }
                         }
                     }
@@ -263,15 +263,15 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
                 }
             }
             sigma[[ri]] <- NA
-            logEC50[[ri]] <- NA
-            stderrlogEC50[[ri]] <- NA
+            logED50[[ri]] <- NA
+            stderrlogED50[[ri]] <- NA
             a[[ri]] <- NA
             b[[ri]] <- NA
             c[[ri]] <- NA
         }
     }
-    results <- data.frame(rsubstance, rn, rlld, rlhd, mtype, logEC50, stderrlogEC50, runit, sigma, a, b)
-    names(results) <- c("Substance","n","lld","lhd","mtype","logEC50","std","unit","sigma","a","b")
+    results <- data.frame(rsubstance, rn, rlld, rlhd, mtype, logED50, stderrlogED50, runit, sigma, a, b)
+    names(results) <- c("Substance","n","lld","lhd","mtype","logED50","std","unit","sigma","a","b")
     if (linlogit) {
         results$c <- c
     }
@@ -303,8 +303,8 @@ drplot <- function(drresults, data, dtype = "std", alpha = 0.95,
         hr <- max(data$response)
         dsubstances <- levels(data$substance)    
     } else {
-        lld <- min(drresults[["logEC50"]],na.rm=TRUE) - 2
-        lhd <- max(drresults[["logEC50"]],na.rm=TRUE) + 2
+        lld <- min(drresults[["logED50"]],na.rm=TRUE) - 2
+        lhd <- max(drresults[["logED50"]],na.rm=TRUE) + 2
         if (length(subset(drresults,mtype=="linlogit")$Substance) != 0) {
             hr <- 1.8 
         } else {
@@ -419,15 +419,15 @@ drplot <- function(drresults, data, dtype = "std", alpha = 0.95,
                 if (nf > 0) {
                     for (j in 1:nf)
                     {
-                        logEC50 <- fits[j,"logEC50"]
+                        logED50 <- fits[j,"logED50"]
                         mtype <- as.character(fits[j, "mtype"])
                         if (mtype == "probit") {
                             scale <- fits[j,"b"]
-                            plot(function(x) pnorm(-x,-logEC50,scale),lld - 0.5, lhd + 2, add=TRUE,col=color)
+                            plot(function(x) pnorm(-x,-logED50,scale),lld - 0.5, lhd + 2, add=TRUE,col=color)
                         }
                         if (mtype == "logit") {
                             scale <- fits[j,"b"]
-                            plot(function(x) plogis(-x,-logEC50,scale),lld - 0.5, lhd + 2, add=TRUE,col=color)
+                            plot(function(x) plogis(-x,-logED50,scale),lld - 0.5, lhd + 2, add=TRUE,col=color)
                         }
                         if (mtype == "weibull") {
                             location <- fits[j,"a"]
@@ -435,7 +435,7 @@ drplot <- function(drresults, data, dtype = "std", alpha = 0.95,
                             plot(function(x) pweibull(-x+location,shape),lld - 0.5, lhd + 2, add=TRUE,col=color)
                         }
                         if (mtype == "linlogit") {
-                            plot(function(x) linlogitf(10^x,1,fits[j,"c"],fits[j,"logEC50"],fits[j,"b"]),
+                            plot(function(x) linlogitf(10^x,1,fits[j,"c"],fits[j,"logED50"],fits[j,"b"]),
                                 lld - 0.5, lhd + 2,
                                 add=TRUE,col=color)
                         }