aboutsummaryrefslogtreecommitdiff
path: root/R/drfit.R
diff options
context:
space:
mode:
authorranke <ranke@d1b72e20-2ee0-0310-a1c4-ad5adbbefcdc>2005-07-16 09:12:22 +0000
committerranke <ranke@d1b72e20-2ee0-0310-a1c4-ad5adbbefcdc>2005-07-16 09:12:22 +0000
commitb83cc0e0ec978b5d047924a222a387152b086f4b (patch)
tree1dff43e042228923c09ef8eea1a078c47f41e26f /R/drfit.R
parentf78a6461cb0510d606fa420bd8fc1683a8815d4b (diff)
Several changes in this revision:
- Addition of the weibull fit - Renaming of lognorm to probit, logis to logit, and linearlogis to linlogit - Fix of a bug that prevented the logis fit to work in recent version - Different way of listing parameters in the output. - Documentation of changes git-svn-id: http://kriemhild.uft.uni-bremen.de/svn/drfit@29 d1b72e20-2ee0-0310-a1c4-ad5adbbefcdc
Diffstat (limited to 'R/drfit.R')
-rw-r--r--R/drfit.R135
1 files changed, 87 insertions, 48 deletions
diff --git a/R/drfit.R b/R/drfit.R
index 38f8e75..603cf94 100644
--- a/R/drfit.R
+++ b/R/drfit.R
@@ -31,14 +31,14 @@ drdata <- function(substances, experimentator = "%", db = "cytotox",
return(data)
}
-linearlogisf <- function(x,k,f,mu,b)
+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,
- lognorm = TRUE, logis = FALSE,
- linearlogis = FALSE, linearlogisWrong = NA, allWrong = NA,
+ probit = TRUE, logit = FALSE, weibull = FALSE,
+ linlogit = FALSE, linlogitWrong = NA, allWrong = NA,
b0 = 2, f0 = 0)
{
if(!is.null(data$ok)) data <- subset(data,ok!="no fit") # Don't use data where ok
@@ -56,9 +56,7 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
sigma <- array() # the standard deviation of the residuals
logEC50 <- vector()
stderrlogEC50 <- vector()
- slope <- vector()
- b <- vector()
- f <- vector()
+ a <- b <- c <- vector()
splitted <- split(data,data$substance)
for (i in substances) {
@@ -95,10 +93,10 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
if (responseathighestdose < 0.5) {
inactive <- FALSE
- if (linearlogis &&
- length(subset(linearlogisWrong,linearlogisWrong == i))==0 &&
+ if (linlogit &&
+ length(subset(linlogitWrong,linlogitWrong == i))==0 &&
length(subset(allWrong,allWrong == i))==0) {
- m <- try(nls(response ~ linearlogisf(dose,1,f,logEC50,b),
+ m <- try(nls(response ~ linlogitf(dose,1,f,logEC50,b),
data=tmp,
start=list(f=f0,logEC50=startlogEC50[[i]],b=b0)))
if (!inherits(m, "try-error")) {
@@ -111,27 +109,28 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
runit[[ri]] <- unit
rlld[[ri]] <- log10(lowestdose)
rlhd[[ri]] <- log10(highestdose)
- mtype[[ri]] <- "linearlogis"
+ mtype[[ri]] <- "linlogit"
logEC50[[ri]] <- coef(m)[["logEC50"]]
- slope[[ri]] <- NA
if (logEC50[[ri]] > rlhd[[ri]]) {
logEC50[[ri]] <- NA
stderrlogEC50[[ri]] <- NA
+ a[[ri]] <- NA
b[[ri]] <- NA
- f[[ri]] <- NA
+ c[[ri]] <- NA
} else {
stderrlogEC50[[ri]] <- s$parameters["logEC50","Std. Error"]
+ a[[ri]] <- coef(m)[["logEC50"]]
b[[ri]] <- coef(m)[["b"]]
- f[[ri]] <- coef(m)[["f"]]
+ c[[ri]] <- coef(m)[["f"]]
}
}
}
- if (logis &&
+ if (weibull &&
length(subset(allWrong,allWrong == i))==0) {
- m <- try(nls(response ~ plogis(-log10(dose),-logEC50,slope),
+ m <- try(nls(response ~ pweibull(-log10(dose)+location,shape),
data=tmp,
- start=list(logEC50=startlogEC50[[i]],slope=1)))
+ start=list(location=-0.16,shape=0.4659)))
if (chooseone==FALSE || fit==FALSE) {
if (!inherits(m, "try-error")) {
fit <- TRUE
@@ -140,29 +139,67 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
sigma[[ri]] <- s$sigma
rsubstance[[ri]] <- i
rn[[ri]] <- n
+ runit[[ri]] <- unit
rlld[[ri]] <- log10(lowestdose)
rlhd[[ri]] <- log10(highestdose)
- mtype[[ri]] <- "logis"
+ mtype[[ri]] <- "weibull"
+ a[[ri]] <- coef(m)[["location"]]
+ b[[ri]] <- coef(m)[["shape"]]
+ sqrdev <- function(logdose) {
+ (0.5 - pweibull( - logdose + a[[ri]], b[[ri]]))^2
+ }
+ logEC50[[ri]] <- nlm(sqrdev,startlogEC50[[i]])$estimate
+ c[[ri]] <- NA
+ if (logEC50[[ri]] > rlhd[[ri]]) {
+ logEC50[[ri]] <- NA
+ stderrlogEC50[[ri]] <- NA
+ a[[ri]] <- NA
+ b[[ri]] <- NA
+ } else {
+ stderrlogEC50[[ri]] <- NA
+ }
+ }
+ }
+ }
+
+ if (logit &&
+ length(subset(allWrong,allWrong == i))==0) {
+ m <- try(nls(response ~ plogis(-log10(dose),-logEC50,scale),
+ data=tmp,
+ start=list(logEC50=startlogEC50[[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
+ rn[[ri]] <- n
+ runit[[ri]] <- unit
+ rlld[[ri]] <- log10(lowestdose)
+ rlhd[[ri]] <- log10(highestdose)
+ mtype[[ri]] <- "logit"
logEC50[[ri]] <- coef(m)[["logEC50"]]
- b[[ri]] <- NA
- f[[ri]] <- NA
+ c[[ri]] <- NA
if (logEC50[[ri]] > rlhd[[ri]]) {
logEC50[[ri]] <- NA
- slope[[ri]] <- NA
stderrlogEC50[[ri]] <- NA
+ a[[ri]] <- NA
+ b[[ri]] <- NA
} else {
- slope[[ri]] <- coef(m)[["slope"]]
stderrlogEC50[[ri]] <- s$parameters["logEC50","Std. Error"]
+ a[[ri]] <- coef(m)[["logEC50"]]
+ b[[ri]] <- coef(m)[["scale"]]
}
}
}
}
- if (lognorm &&
+ if (probit &&
length(subset(allWrong,allWrong == i))==0) {
- m <- try(nls(response ~ pnorm(-log10(dose),-logEC50,slope),
+ m <- try(nls(response ~ pnorm(-log10(dose),-logEC50,scale),
data=tmp,
- start=list(logEC50=startlogEC50[[i]],slope=1)))
+ start=list(logEC50=startlogEC50[[i]],scale=1)))
if (chooseone==FALSE || fit==FALSE) {
if (!inherits(m, "try-error")) {
fit <- TRUE
@@ -174,17 +211,18 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
runit[[ri]] <- unit
rlld[[ri]] <- log10(lowestdose)
rlhd[[ri]] <- log10(highestdose)
- mtype[[ri]] <- "lognorm"
+ mtype[[ri]] <- "probit"
logEC50[[ri]] <- coef(m)[["logEC50"]]
- b[[ri]] <- NA
- f[[ri]] <- NA
+ c[[ri]] <- NA
if (logEC50[[ri]] > rlhd[[ri]]) {
logEC50[[ri]] <- NA
- slope[[ri]] <- NA
stderrlogEC50[[ri]] <- NA
+ a[[ri]] <- NA
+ b[[ri]] <- NA
} else {
- slope[[ri]] <- coef(m)[["slope"]]
stderrlogEC50[[ri]] <- s$parameters["logEC50","Std. Error"]
+ a[[ri]] <- coef(m)[["logEC50"]]
+ b[[ri]] <- coef(m)[["scale"]]
}
}
}
@@ -215,19 +253,15 @@ drfit <- function(data, startlogEC50 = NA, chooseone=TRUE,
sigma[[ri]] <- NA
logEC50[[ri]] <- NA
stderrlogEC50[[ri]] <- NA
- slope[[ri]] <- NA
+ a[[ri]] <- NA
b[[ri]] <- NA
- f[[ri]] <- NA
+ c[[ri]] <- NA
}
}
- results <- data.frame(rsubstance, rn, rlld, rlhd, mtype, logEC50, stderrlogEC50, runit, sigma)
- names(results) <- c("Substance","n","lld","lhd","mtype","logEC50","std","unit","sigma")
- if (lognorm || logis) {
- results$slope <- slope
- }
- if (linearlogis) {
- results$b <- b
- results$f <- f
+ 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")
+ if (linlogit) {
+ results$c <- c
}
return(results)
}
@@ -257,7 +291,7 @@ drplot <- function(drresults, data, dtype = "std", alpha = 0.95,
} else {
lld <- min(drresults[["logEC50"]],na.rm=TRUE) - 2
lhd <- max(drresults[["logEC50"]],na.rm=TRUE) + 2
- if (length(subset(drresults,mtype=="linearlogis")$Substance) != 0) {
+ if (length(subset(drresults,mtype=="linlogit")$Substance) != 0) {
hr <- 1.8
} else {
hr <- 1.0
@@ -372,16 +406,21 @@ drplot <- function(drresults, data, dtype = "std", alpha = 0.95,
{
logEC50 <- fits[j,"logEC50"]
mtype <- as.character(fits[j, "mtype"])
- if (mtype == "lognorm") {
- slope <- fits[j,"slope"]
- plot(function(x) pnorm(-x,-logEC50,slope),lld - 0.5, lhd + 2, add=TRUE,col=color)
+ if (mtype == "probit") {
+ scale <- fits[j,"b"]
+ plot(function(x) pnorm(-x,-logEC50,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)
}
- if (mtype == "logis") {
- slope <- fits[j,"slope"]
- plot(function(x) plogis(-x,-logEC50,slope),lld - 0.5, lhd + 2, add=TRUE,col=color)
+ if (mtype == "weibull") {
+ location <- fits[j,"a"]
+ shape <- fits[j,"b"]
+ plot(function(x) pweibull(-x+location,shape),lld - 0.5, lhd + 2, add=TRUE,col=color)
}
- if (mtype == "linearlogis") {
- plot(function(x) linearlogisf(10^x,1,fits[j,"f"],fits[j,"logEC50"],fits[j,"b"]),
+ if (mtype == "linlogit") {
+ plot(function(x) linlogitf(10^x,1,fits[j,"c"],fits[j,"logEC50"],fits[j,"b"]),
lld - 0.5, lhd + 2,
add=TRUE,col=color)
}

Contact - Imprint