Cleaned up version, only containing very basic stuff.

git-svn-id: http://kriemhild.uft.uni-bremen.de/svn/chemCal@2 5fad18fb-23f0-0310-ab10-e59a3bee62b4

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)
 }