A trial to improve the lod function, I don't really understand

the predict.lm function unfortunately



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

6 files changed, 26 insertions, 21 deletions

diff --git a/R/inverse.predict.lm.R b/R/inverse.predict.lm.R
index b48c967..2ead9fb 100644
--- a/R/inverse.predict.lm.R
+++ b/R/inverse.predict.lm.R
@@ -55,11 +55,12 @@ inverse.predict.rlm <- function(object, newdata, ...,
 
   yx <- split(object$model$y,object$model$x)
   n <- length(yx)
+  df <- n - length(objects$coef)
   x <- as.numeric(names(yx))
   ybar <- sapply(yx,mean)
   yhatx <- split(object$fitted.values,object$model$x)
   yhat <- sapply(yhatx,mean)
-  se <- sqrt(sum(w*(ybar - yhat)^2)/(n-2))
+  se <- sqrt(sum(w*(ybar - yhat)^2)/df)
   if (ss == "auto") {
     ss <- se
   } else {
diff --git a/R/lod.R b/R/lod.R
index 6562c04..1bb7981 100644
--- a/R/lod.R
+++ b/R/lod.R
@@ -1,23 +1,27 @@
-lod <- function(object, ..., alpha = 0.05, beta = 0.05, n = 1, w = "auto")
+lod <- function(object, ..., alpha = 0.05, beta = 0.05)
 {
   UseMethod("lod")
 }
 
-lod.default <- function(object, ..., alpha = 0.05, beta = 0.05,
-  n = 1, w = "auto")
+lod.default <- function(object, ..., alpha = 0.05, beta = 0.05)
 {
   stop("lod is only implemented for univariate lm objects.")
 }
 
-lod.lm <- function(object, ..., alpha = 0.05, beta = 0.05, n = 1, w = "auto")
+lod.lm <- function(object, ..., alpha = 0.05, beta = 0.05)
 {
-  f <- function(x) {
-    y1 <- predict(object, data.frame(x = 0), interval="prediction", 
+  y0 <- predict(object, data.frame(x = 0), interval="prediction", 
       level = 1 - alpha)
-    y2 <- predict(object, data.frame(x = x), interval="prediction",
-      level = 1 - beta)
-    (y2[[1,"lwr"]] - y1[[1,"upr"]])^2
+  yc <- y0[[1,"upr"]]
+  xc <- inverse.predict(object,yc)[["Prediction"]]
+  f <- function(x)
+  {
+      # Here I need the variance of y values as a function of x or
+      # y values
+      # Strangely, setting the confidence level to 0.5 does not result
+      # in a zero confidence or prediction interval
   }
-  tmp <- optimize(f,interval=c(0,max(object$model$x)))
-  return(tmp$minimum)
+  lod.x <- optimize(f,interval=c(0,max(object$model$x)))$minimum
+  lod.y <-  predict(object, data.frame(x = lod.x))
+  return(list(x = lod.x, y = lod.y))
 }
diff --git a/inst/doc/chemCal-001.pdf b/inst/doc/chemCal-001.pdf
index e2992fc..8375fd5 100644
--- a/inst/doc/chemCal-001.pdf
+++ b/inst/doc/chemCal-001.pdf
@@ -2,8 +2,8 @@
 %���ρ�\r
 1 0 obj
 <<
-/CreationDate (D:20060512235955)
-/ModDate (D:20060512235955)
+/CreationDate (D:20060515234125)
+/ModDate (D:20060515234125)
 /Title (R Graphics Output)
 /Producer (R 2.3.0)
 /Creator (R)
diff --git a/inst/doc/chemCal.log b/inst/doc/chemCal.log
index 110dbec..785a8cf 100644
--- a/inst/doc/chemCal.log
+++ b/inst/doc/chemCal.log
@@ -1,4 +1,4 @@
-This is pdfeTeX, Version 3.141592-1.21a-2.2 (Web2C 7.5.4) (format=pdflatex 2006.4.5)  12 MAY 2006 23:59
+This is pdfeTeX, Version 3.141592-1.21a-2.2 (Web2C 7.5.4) (format=pdflatex 2006.4.8)  15 MAY 2006 23:41
 entering extended mode
 **chemCal.tex
 (./chemCal.tex
@@ -325,8 +325,8 @@ LaTeX Font Info:    External font `cmex10' loaded for size
  [2 <./chemCal-001.pdf>]
 [3] (./chemCal.aux) ) 
 Here is how much of TeX's memory you used:
- 3764 strings out of 94499
- 51832 string characters out of 1175814
+ 3764 strings out of 94500
+ 51832 string characters out of 1175772
  96290 words of memory out of 1000000
  6863 multiletter control sequences out of 10000+50000
  21774 words of font info for 52 fonts, out of 500000 for 2000
@@ -346,4 +346,4 @@ type1/bluesky/cm/cmbx12.pfb> </var/cache/fonts/pk/ljfour/jknappen/tc/tctt1000.6
 xmf-tetex/fonts/type1/bluesky/cm/cmtt10.pfb></usr/share/texmf-tetex/fonts/type1
 /bluesky/cm/cmr10.pfb></usr/share/texmf-tetex/fonts/type1/bluesky/cm/cmr12.pfb>
 </usr/share/texmf-tetex/fonts/type1/bluesky/cm/cmr17.pfb>
-Output written on chemCal.pdf (3 pages, 107518 bytes).
+Output written on chemCal.pdf (3 pages, 107539 bytes).
diff --git a/inst/doc/chemCal.pdf b/inst/doc/chemCal.pdf
index 8e3249b..e73cc52 100644
--- a/inst/doc/chemCal.pdf
+++ b/inst/doc/chemCal.pdf
diff --git a/inst/doc/chemCal.tex b/inst/doc/chemCal.tex
index 4dbdce3..a9e8495 100644
--- a/inst/doc/chemCal.tex
+++ b/inst/doc/chemCal.tex
@@ -63,13 +63,13 @@ $Prediction
 [1] 5.865367
 
 $`Standard Error`
-[1] 0.892611
+[1] 0.7288138
 
 $Confidence
-[1] 2.478285
+[1] 2.023511
 
 $`Confidence Limits`
-[1] 3.387082 8.343652
+[1] 3.841856 7.888878
 \end{Soutput}
 \end{Schunk}