aboutsummaryrefslogtreecommitdiff
path: root/R/chemCal.R
diff options
context:
space:
mode:
Diffstat (limited to 'R/chemCal.R')
-rw-r--r--R/chemCal.R132
1 files changed, 16 insertions, 116 deletions
diff --git a/R/chemCal.R b/R/chemCal.R
index 5de33f8..fab7db4 100644
--- a/R/chemCal.R
+++ b/R/chemCal.R
@@ -20,10 +20,11 @@ calm <- function(data)
predict.calm <- predict.lm
print.calm <- print.lm
summary.calm <- summary.lm
-plot.calm <- function(m,
+plot.calm <- function(x,...,
xunit="",yunit="",measurand="",
level=0.95)
{
+ m <- x
x <- m$model$x
y <- m$model$y
newdata <- data.frame(x = seq(0,max(x),length=250))
@@ -59,137 +60,36 @@ plot.calm <- function(m,
col = c("black", "green4", "red"),
horiz = FALSE, cex = 0.9, bg = "gray95")
}
-predictx.calm <- function(m,measurements)
+predictx <- function(m,yobs,level=0.95)
{
s <- summary(m)
xi <- m$model$x
yi <- m$model$y
n <- length(yi)
- yobs <- newdata[[1]]
p <- length(yobs)
- if (!m$intercept)
- {
- varb1 <- summary(m)$coef["x","Std. Error"]^2
- xpred <- mean(yobs)/b1
- varxpred <- (varyobs + xpred^2 * varb1) / b1^2
- sdxpred <- sqrt(varxpred)
- } else
- {
- b0 <- summary(m)$coef["(Intercept)","Estimate"]
- b1 <- summary(m)$coef["xi","Estimate"]
- xpred <- (mean(yobs) - b0)/b1
- sumxxbar <- sum((xi - mean(xi))^2)
- if (!syobs)
- {
- yybar <- (mean(yobs) - mean(yi))^2
- sdxpred <- (S/b1) * (1/p + 1/n + yybar/(b1^2 * sumxxbar))^0.5
- } else
- {
- sdxpred <- ((varyobs^0.5/b1) + (S/b1)^2 * (1/n + ((xpred - mean(xi))^2)/sumxxbar))^0.5
- }
- }
- t <- qt((1 + level)/2,ntot - 2)
- confxpred <- t * sdxpred
-
- result <- c(estimate=xpred,sdxpred=sdxpred,syobs=syobs,
- confxpred=confxpred)
- digits <- max(c(3,round(log10(xpred/confxpred)) + 2))
- roundedresult <- round(result,digits=digits)
- confidenceinterval <- paste(roundedresult["estimate"],"+-",
- roundedresult["confxpred"],xunit)
- roundedresult[["confidenceinterval"]] <- confidenceinterval
- invisible(roundedresult)
-}
-calpredict <- function(yobs,xi,yi,xunit="",level=0.95,intercept=FALSE,syobs=FALSE)
-{
- if (length(xi)!=length(yi))
- {
- cat("xi and yi have to be of the same length\n")
- }
- xi <- xi[!is.na(yi)]
- yi <- yi[!is.na(yi)]
- n <- length(yi)
- p <- length(yobs)
- if (!intercept)
- {
- m <- lm(yi ~ xi - 1)
- } else
- {
- m <- lm(yi ~ xi)
- }
- S <- summary(m)$sigma
- b1 <- summary(m)$coef["xi","Estimate"]
-
- if (syobs)
- {
- if (is.numeric(syobs))
- {
- varyobs <- syobs^2
- ntot <- n
- } else
- {
- if (length(yobs)==1)
- {
- cat("yobs has to contain more than one number vector, if you use syobs=TRUE\n")
- }
- varyobs <- var(yobs)
- ntot <- n + p
- }
- } else
- {
- varyobs <- S^2
- ntot <- n
+ if (p > 1) {
+ varyobs <- var(yobs)
+ } else {
+ varyobs <- 0
}
-
- if (!intercept)
- {
- varb1 <- summary(m)$coef["xi","Std. Error"]^2
+ if (!m$intercept) {
+ b1 <- summary(m)$coef["x","Estimate"]
+ varb1 <- summary(m)$coef["x","Std. Error"]^2
xpred <- mean(yobs)/b1
varxpred <- (varyobs + xpred^2 * varb1) / b1^2
sdxpred <- sqrt(varxpred)
} else
{
b0 <- summary(m)$coef["(Intercept)","Estimate"]
- b1 <- summary(m)$coef["xi","Estimate"]
+ b1 <- summary(m)$coef["x","Estimate"]
+ S <- summary(m)$sigma
xpred <- (mean(yobs) - b0)/b1
sumxxbar <- sum((xi - mean(xi))^2)
- if (!syobs)
- {
- yybar <- (mean(yobs) - mean(yi))^2
- sdxpred <- (S/b1) * (1/p + 1/n + yybar/(b1^2 * sumxxbar))^0.5
- } else
- {
- sdxpred <- ((varyobs^0.5/b1) + (S/b1)^2 * (1/n + ((xpred - mean(xi))^2)/sumxxbar))^0.5
- }
+ yybar <- (mean(yobs) - mean(yi))^2
+ sdxpred <- (S/b1) * (1/p + 1/n + yybar/(b1^2 * sumxxbar))^0.5
}
- t <- qt((1 + level)/2,ntot - 2)
+ t <- qt((1 + level)/2,n - 2)
confxpred <- t * sdxpred
- result <- c(estimate=xpred,sdxpred=sdxpred,syobs=syobs,
- confxpred=confxpred)
- digits <- max(c(3,round(log10(xpred/confxpred)) + 2))
- roundedresult <- round(result,digits=digits)
- confidenceinterval <- paste(roundedresult["estimate"],"+-",
- roundedresult["confxpred"],xunit)
- roundedresult[["confidenceinterval"]] <- confidenceinterval
- invisible(roundedresult)
-}
-
-multical <- function(cf,df,intercept=FALSE)
-{
- rf <- data.frame(name=levels(df$name))
- substances <- colnames(df)[-1]
- for (s in substances)
- {
- r <- vector()
- for (sample in levels(df$name))
- {
- tmp <- calpredict(subset(df,name==sample)[[s]],
- cf[["conc"]],cf[[s]],
- intercept=intercept)
- r <- c(r,tmp[["confidenceinterval"]])
- }
- rf[[s]] <- r
- }
- return(rf)
+ result <- c(estimate=xpred,sdxpred=sdxpred,confxpred=confxpred)
}

Contact - Imprint