aboutsummaryrefslogtreecommitdiff
path: root/R
diff options
context:
space:
mode:
Diffstat (limited to 'R')
-rw-r--r--R/drfit.R26
-rw-r--r--R/drplot.R30
2 files changed, 41 insertions, 15 deletions
diff --git a/R/drfit.R b/R/drfit.R
index ff03d90..1bbffa8 100644
--- a/R/drfit.R
+++ b/R/drfit.R
@@ -1,6 +1,7 @@
drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
probit = TRUE, logit = FALSE, weibull = FALSE,
- linlogit = FALSE, linlogitWrong = NA, allWrong = NA,
+ linlogit = FALSE, conf = FALSE,
+ 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
@@ -20,6 +21,7 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
sigma <- array() # the standard deviation of the residuals
logED50 <- vector()
stderrlogED50 <- vector()
+ conflogED50 <- vector()
a <- b <- c <- vector()
splitted <- split(data,data$substance)
@@ -47,7 +49,7 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
n <- ndl <- 0
} else {
ndl <- length(levels(factor(tmp$dose)))
- n <- round(length(tmp$response)/ndl)
+ n <- length(tmp$response)
if (is.na(startlogED50[i])){
w <- 1/abs(tmp$response - 0.3)
startlogED50[[i]] <- sum(w * log10(tmp$dose))/sum(w)
@@ -86,12 +88,14 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
mtype[[ri]] <- "no fit"
logED50[[ri]] <- NA
stderrlogED50[[ri]] <- NA
+ conflogED50[[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)
a[[ri]] <- coef(m)[["logED50"]]
b[[ri]] <- coef(m)[["b"]]
c[[ri]] <- coef(m)[["f"]]
@@ -122,11 +126,13 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
mtype[[ri]] <- "no fit"
logED50[[ri]] <- NA
stderrlogED50[[ri]] <- NA
+ conflogED50[[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)
a[[ri]] <- coef(m)[["logED50"]]
b[[ri]] <- coef(m)[["scale"]]
}
@@ -158,11 +164,13 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
mtype[[ri]] <- "no fit"
logED50[[ri]] <- NA
stderrlogED50[[ri]] <- NA
+ conflogED50[[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)
}
}
}
@@ -196,11 +204,13 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
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)
}
}
}
@@ -238,13 +248,21 @@ drfit <- function(data, startlogED50 = NA, chooseone=TRUE,
sigma[[ri]] <- NA
logED50[[ri]] <- NA
stderrlogED50[[ri]] <- NA
+ conflogED50[[ri]] <- NA
a[[ri]] <- NA
b[[ri]] <- NA
c[[ri]] <- NA
}
}
- 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")
+ 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")
+ }
+
if (linlogit) {
results$c <- c
}
diff --git a/R/drplot.R b/R/drplot.R
index baec69d..8faaa0f 100644
--- a/R/drplot.R
+++ b/R/drplot.R
@@ -1,6 +1,7 @@
drplot <- function(drresults, data,
dtype = "std", alpha = 0.95, ctype = "none",
path = "./", fileprefix = "drplot", overlay = FALSE,
+ xlim = c("auto","auto"), ylim = c("auto","auto"),
postscript = FALSE, pdf = FALSE, png = FALSE,
bw = TRUE,
pointsize = 12,
@@ -14,14 +15,15 @@ drplot <- function(drresults, data,
unit <- "different units"
}
- # Get the plot limits on the x-axis (log of the dose) and y axis
+ # Determine the plot limits on the x-axis and y axis
if(is.data.frame(data)) {
+ # Get rid of pseudo substance names of controls
nonzerodata <- subset(data,dose!=0)
- nonzerodata$substance <- factor(nonzerodata$substance) # Get rid of pseudo substance names of controls
+ nonzerodata$substance <- factor(nonzerodata$substance)
zerodata <- subset(data,dose==0)
nc <- length(zerodata$dose) # Number of control points
if (nc > 0) {
- sdc <- sd(zerodata$response) # Standard deviation of control responses
+ sdc <- sd(zerodata$response)
controlconf <- sdc * qt((1 + alpha)/2, nc - 1) / sqrt(nc)
cat("There are ",nc,"data points with dose 0 (control values)\n")
cat("with a standard deviation of",sdc,"\n")
@@ -52,6 +54,12 @@ drplot <- function(drresults, data,
hr <- 1.0
}
}
+ if (xlim[1] == "auto") xlim[1] <- lld - 0.5
+ if (xlim[2] == "auto") xlim[2] <- lhd + 1
+ if (ylim[1] == "auto") ylim[1] <- -0.1
+ if (ylim[2] == "auto") ylim[2] <- hr + 0.2
+ xlim <- as.numeric(xlim)
+ ylim <- as.numeric(ylim)
# Prepare overlay plot if requested
if (overlay)
@@ -79,10 +87,10 @@ drplot <- function(drresults, data,
}
plot(0,type="n",
- xlim=c(lld - 0.5, lhd + 1),
- ylim= c(-0.1, hr + 0.2),
- xlab=paste("Decadic Logarithm of the dose in ", unit),
- ylab="Normalized response")
+ xlim = xlim,
+ ylim = ylim,
+ xlab = paste("Decadic Logarithm of the dose in ", unit),
+ ylab = "Normalized response")
}
# Plot the data either as raw data or as error bars
@@ -123,10 +131,10 @@ drplot <- function(drresults, data,
}
plot(0,type="n",
- xlim=c(lld - 0.5, lhd + 2),
- ylim= c(-0.1, hr + 0.2),
- xlab=paste("Decadic Logarithm of the dose in ", unit),
- ylab="Normalized response")
+ xlim = xlim,
+ ylim = ylim,
+ xlab = paste("Decadic Logarithm of the dose in ", unit),
+ ylab = "Normalized response")
}
if (!overlay) legend(lpos, i,lty = 1, col = color, inset=0.05)
tmp$dosefactor <- factor(tmp$dose) # necessary because the old

Contact - Imprint