From 424d3633234125cf2c8f65fd220fb92ddc0d2686 Mon Sep 17 00:00:00 2001 From: ranke Date: Mon, 27 Feb 2006 10:49:58 +0000 Subject: 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 --- R/drfit.R | 227 ++++++++++++++++++++++++++++++++------------------------------ 1 file changed, 118 insertions(+), 109 deletions(-) (limited to 'R') 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 -- cgit v1.2.1