Now the model type (mtype) is set to "active" if the response at the lowest

dose level is smaller than 0.5.



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

1 files changed, 118 insertions, 109 deletions

diff --git a/R/drfit.R b/R/drfit.R
index 9eee87a..1488cd4 100644
--- a/R/drfit.R
+++ b/R/drfit.R
@@ -101,50 +101,19 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
             lhd <- log10(highestdose)
             lld <- log10(lowestdose)
             responseathighestdose <- mean(subset(tmp,dose==highestdose)$response)
+            responseatlowestdose <- mean(subset(tmp,dose==lowestdose)$response)
             if (responseathighestdose < 0.5) {
                 inactive <- FALSE
-
-                if (linlogit && 
-                    length(subset(linlogitWrong,linlogitWrong == i))==0 &&
-                    length(subset(allWrong,allWrong == i))==0) {
-                    m <- try(nls(response ~ linlogitf(dose,1,f,logED50,b),
-                            data=tmp,
-                            start=list(f=f0,logED50=startlogED50[[i]],b=b0)))
-                    if (!inherits(m, "try-error")) {
-                        fit <- TRUE
-                        ri <- ri + 1
-                        s <- summary(m)
-                        sigma[[ri]] <- s$sigma
-                        rsubstance[[ri]] <- i
-                        rndl[[ri]] <- ndl
-                        rn[[ri]] <- n
-                        runit[[ri]] <- unit
-                        rlld[[ri]] <- log10(lowestdose)
-                        rlhd[[ri]] <- log10(highestdose)
-                        logED50[[ri]] <- coef(m)[["logED50"]]
-                        if (logED50[[ri]] > rlhd[[ri]]) {
-                            mtype[[ri]] <- "no fit"
-                            logED50[[ri]] <- NA
-                            stderrlogED50[[ri]] <- NA
-                            a[[ri]] <- NA
-                            b[[ri]] <- NA
-                            c[[ri]] <- NA
-                        } else {
-                            mtype[[ri]] <- "linlogit"
-                            stderrlogED50[[ri]] <- s$parameters["logED50","Std. Error"]
-                            a[[ri]] <- coef(m)[["logED50"]]
-                            b[[ri]] <- coef(m)[["b"]]
-                            c[[ri]] <- coef(m)[["f"]]
-                        }
-                    }
-                }
-
-                if (probit &&
-                    length(subset(allWrong,allWrong == i))==0) {
-                    m <- try(nls(response ~ pnorm(-log10(dose),-logED50,scale),
+                if (responseatlowestdose < 0.5) {
+                    active <- TRUE
+                } else {
+                    active <- FALSE
+                    if (linlogit && 
+                        length(subset(linlogitWrong,linlogitWrong == i))==0 &&
+                        length(subset(allWrong,allWrong == i))==0) {
+                        m <- try(nls(response ~ linlogitf(dose,1,f,logED50,b),
                                 data=tmp,
-                                start=list(logED50=startlogED50[[i]],scale=1)))
-                    if (chooseone==FALSE || fit==FALSE) {
+                                start=list(f=f0,logED50=startlogED50[[i]],b=b0)))
                         if (!inherits(m, "try-error")) {
                             fit <- TRUE
                             ri <- ri + 1
@@ -157,95 +126,131 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
                             rlld[[ri]] <- log10(lowestdose)
                             rlhd[[ri]] <- log10(highestdose)
                             logED50[[ri]] <- coef(m)[["logED50"]]
-                            c[[ri]] <- NA
                             if (logED50[[ri]] > rlhd[[ri]]) {
                                 mtype[[ri]] <- "no fit"
                                 logED50[[ri]] <- NA
                                 stderrlogED50[[ri]] <- NA
                                 a[[ri]] <- NA
                                 b[[ri]] <- NA
+                                c[[ri]] <- NA
                             } else {
-                                mtype[[ri]] <- "probit"
+                                mtype[[ri]] <- "linlogit"
                                 stderrlogED50[[ri]] <- s$parameters["logED50","Std. Error"]
                                 a[[ri]] <- coef(m)[["logED50"]]
-                                b[[ri]] <- coef(m)[["scale"]]
+                                b[[ri]] <- coef(m)[["b"]]
+                                c[[ri]] <- coef(m)[["f"]]
                             }
                         }
                     }
-                }
 
-                if (logit &&
-                    length(subset(allWrong,allWrong == i))==0) {
-                    m <- try(nls(response ~ plogis(-log10(dose),-logED50,scale),
-                            data=tmp,
-                            start=list(logED50=startlogED50[[i]],scale=1)))
-                    if (chooseone==FALSE || fit==FALSE) {
-                        if (!inherits(m, "try-error")) {
-                            fit <- TRUE
-                            ri <- ri + 1
-                            s <- summary(m)
-                            sigma[[ri]] <- s$sigma
-                            rsubstance[[ri]] <- i
-                            rndl[[ri]] <- ndl
-                            rn[[ri]] <- n
-                            runit[[ri]] <- unit
-                            rlld[[ri]] <- log10(lowestdose)
-                            rlhd[[ri]] <- log10(highestdose)
-                            logED50[[ri]] <- a[[ri]] <- coef(m)[["logED50"]]
-                            b[[ri]] <- coef(m)[["scale"]]
-                            c[[ri]] <- NA
-                            if (logED50[[ri]] > rlhd[[ri]]) {
-                                mtype[[ri]] <- "no fit"
-                                logED50[[ri]] <- NA
-                                stderrlogED50[[ri]] <- NA
-                                a[[ri]] <- NA
-                                b[[ri]] <- NA
-                            } else {
-                                mtype[[ri]] <- "logit"
-                                stderrlogED50[[ri]] <- s$parameters["logED50","Std. Error"]
+                    if (probit &&
+                        length(subset(allWrong,allWrong == i))==0) {
+                        m <- try(nls(response ~ pnorm(-log10(dose),-logED50,scale),
+                                    data=tmp,
+                                    start=list(logED50=startlogED50[[i]],scale=1)))
+                        if (chooseone==FALSE || fit==FALSE) {
+                            if (!inherits(m, "try-error")) {
+                                fit <- TRUE
+                                ri <- ri + 1
+                                s <- summary(m)
+                                sigma[[ri]] <- s$sigma
+                                rsubstance[[ri]] <- i
+                                rndl[[ri]] <- ndl
+                                rn[[ri]] <- n
+                                runit[[ri]] <- unit
+                                rlld[[ri]] <- log10(lowestdose)
+                                rlhd[[ri]] <- log10(highestdose)
+                                logED50[[ri]] <- coef(m)[["logED50"]]
+                                c[[ri]] <- NA
+                                if (logED50[[ri]] > rlhd[[ri]]) {
+                                    mtype[[ri]] <- "no fit"
+                                    logED50[[ri]] <- NA
+                                    stderrlogED50[[ri]] <- NA
+                                    a[[ri]] <- NA
+                                    b[[ri]] <- NA
+                                } else {
+                                    mtype[[ri]] <- "probit"
+                                    stderrlogED50[[ri]] <- s$parameters["logED50","Std. Error"]
+                                    a[[ri]] <- coef(m)[["logED50"]]
+                                    b[[ri]] <- coef(m)[["scale"]]
+                                }
                             }
                         }
                     }
-                }
 
-                if (weibull &&
-                    length(subset(allWrong,allWrong == i))==0) {
-                    m <- try(nls(response ~ pweibull(-log10(dose)+location,shape),
-                            data=tmp,
-                            start=list(location=startlogED50[[i]],shape=s0)))
-                    if (chooseone==FALSE || fit==FALSE) {
-                        if (!inherits(m, "try-error")) {
-                            fit <- TRUE
-                            ri <- ri + 1
-                            s <- summary(m)
-                            sigma[[ri]] <- s$sigma
-                            rsubstance[[ri]] <- i
-                            rndl[[ri]] <- ndl
-                            rn[[ri]] <- n
-                            runit[[ri]] <- unit
-                            rlld[[ri]] <- log10(lowestdose)
-                            rlhd[[ri]] <- log10(highestdose)
-                            a[[ri]] <- coef(m)[["location"]]
-                            b[[ri]] <- coef(m)[["shape"]]
-                            sqrdev <- function(logdose) {
-                                (0.5 - pweibull( - logdose + a[[ri]], b[[ri]]))^2 
+                    if (logit &&
+                        length(subset(allWrong,allWrong == i))==0) {
+                        m <- try(nls(response ~ plogis(-log10(dose),-logED50,scale),
+                                data=tmp,
+                                start=list(logED50=startlogED50[[i]],scale=1)))
+                        if (chooseone==FALSE || fit==FALSE) {
+                            if (!inherits(m, "try-error")) {
+                                fit <- TRUE
+                                ri <- ri + 1
+                                s <- summary(m)
+                                sigma[[ri]] <- s$sigma
+                                rsubstance[[ri]] <- i
+                                rndl[[ri]] <- ndl
+                                rn[[ri]] <- n
+                                runit[[ri]] <- unit
+                                rlld[[ri]] <- log10(lowestdose)
+                                rlhd[[ri]] <- log10(highestdose)
+                                logED50[[ri]] <- a[[ri]] <- coef(m)[["logED50"]]
+                                b[[ri]] <- coef(m)[["scale"]]
+                                c[[ri]] <- NA
+                                if (logED50[[ri]] > rlhd[[ri]]) {
+                                    mtype[[ri]] <- "no fit"
+                                    logED50[[ri]] <- NA
+                                    stderrlogED50[[ri]] <- NA
+                                    a[[ri]] <- NA
+                                    b[[ri]] <- NA
+                                } else {
+                                    mtype[[ri]] <- "logit"
+                                    stderrlogED50[[ri]] <- s$parameters["logED50","Std. Error"]
+                                }
                             }
-                            logED50[[ri]] <- nlm(sqrdev,startlogED50[[i]])$estimate
-                            c[[ri]] <- NA
-                            if (logED50[[ri]] > rlhd[[ri]]) {
-                                mtype[[ri]] <- "no fit"
-                                logED50[[ri]] <- NA
-                                stderrlogED50[[ri]] <- NA
-                                a[[ri]] <- NA
-                                b[[ri]] <- NA
-                            } else {
-                                mtype[[ri]] <- "weibull"
-                                stderrlogED50[[ri]] <- NA
+                        }
+                    }
+
+                    if (weibull &&
+                        length(subset(allWrong,allWrong == i))==0) {
+                        m <- try(nls(response ~ pweibull(-log10(dose)+location,shape),
+                                data=tmp,
+                                start=list(location=startlogED50[[i]],shape=s0)))
+                        if (chooseone==FALSE || fit==FALSE) {
+                            if (!inherits(m, "try-error")) {
+                                fit <- TRUE
+                                ri <- ri + 1
+                                s <- summary(m)
+                                sigma[[ri]] <- s$sigma
+                                rsubstance[[ri]] <- i
+                                rndl[[ri]] <- ndl
+                                rn[[ri]] <- n
+                                runit[[ri]] <- unit
+                                rlld[[ri]] <- log10(lowestdose)
+                                rlhd[[ri]] <- log10(highestdose)
+                                a[[ri]] <- coef(m)[["location"]]
+                                b[[ri]] <- coef(m)[["shape"]]
+                                sqrdev <- function(logdose) {
+                                    (0.5 - pweibull( - logdose + a[[ri]], b[[ri]]))^2 
+                                }
+                                logED50[[ri]] <- nlm(sqrdev,startlogED50[[i]])$estimate
+                                c[[ri]] <- NA
+                                if (logED50[[ri]] > rlhd[[ri]]) {
+                                    mtype[[ri]] <- "no fit"
+                                    logED50[[ri]] <- NA
+                                    stderrlogED50[[ri]] <- NA
+                                    a[[ri]] <- NA
+                                    b[[ri]] <- NA
+                                } else {
+                                    mtype[[ri]] <- "weibull"
+                                    stderrlogED50[[ri]] <- NA
+                                }
                             }
                         }
                     }
-                }
 
+                }
 
             } else {
                 inactive <- TRUE
@@ -267,7 +272,11 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
                 if (inactive) {
                     mtype[[ri]] <- "inactive"
                 } else {
-                    mtype[[ri]] <- "no fit"
+                    if (active) {
+                        mtype[[ri]] <- "active"
+                    } else {
+                        mtype[[ri]] <- "no fit"
+                    }
                 }
             }
             sigma[[ri]] <- NA