diff options
Diffstat (limited to 'R/chemCal.R')
| -rw-r--r-- | R/chemCal.R | 132 |
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) } |