From 0973370a6e27952df81aaae2a05104195e3bf633 Mon Sep 17 00:00:00 2001
From: ranke <ranke@5fad18fb-23f0-0310-ab10-e59a3bee62b4>
Date: Mon, 15 May 2006 23:07:31 +0000
Subject: 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
---
 R/inverse.predict.lm.R |  3 ++-
 R/lod.R                | 26 +++++++++++++++-----------
 2 files changed, 17 insertions(+), 12 deletions(-)

(limited to 'R')

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))
 }
-- 
cgit v1.2.1