From: Ted Harding <Ted.Harding_at_nessie.mcc.ac.uk>

Date: Wed 20 Apr 2005 - 17:58:37 EST

E-Mail: (Ted Harding) <Ted.Harding@nessie.mcc.ac.uk> Fax-to-email: +44 (0)870 094 0861

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Wed Apr 20 18:58:32 2005

Date: Wed 20 Apr 2005 - 17:58:37 EST

Sorry, I was doing this too late last night!
All stands as before except for the calculation at the end
which is now corrected as follows:

On 19-Apr-05 Ted Harding wrote:

[code repeated for reference]

> The following function implements the above expressions.

*> It is a very crude approach to solving the problem, and
**> I'm sure that a more thoughtful approach would lead more
**> directly to the answer, but at least it gets you there
**> eventually.
**>
**> ===========================================
**>
**> R.calib <- function(x,y,X,Y){
**> n<-length(x) ; mx<-mean(x) ; my<-mean(y) ;
**> x<-(x-mx) ; y<-(y-my) ; X<-(X-mx) ; Y<-(Y-my)
**>
**> ah<-0 ; bh<-(sum(x*y))/(sum(x^2)) ; Xh <- Y/bh
**> sh2 <- (sum((y-ah-bh*x)^2))/(n+1)
**>
**> D<-(n+1)*sum(x^2) + n*X^2
**> at<-(Y*sum(x^2) - X*sum(x*y))/D; bt<-((n+1)*sum(x*y) + n*X*Y)/D
**> st2<-(sum((y - at - bt*x)^2) + (Y - at - bt*X)^2)/(n+1)
**>
**> R<-(sh2/st2)
**>
**> F<-(n-2)*(1-R)/R
**>
**> x<-(x+mx) ; y<-(y+my) ;
**> X<-(X+mx) ; Y<-(Y+my) ; Xh<-(Xh+mx) ;
**> PF<-(pf(F,1,(n-2)))
**> list(x=x,y=y,X=X,Y=Y,R=R,F=F,PF=PF,
**> ahat=ah,bhat=bh,sh2=sh2,
**> atil=at,btil=bt,st2=st2,
**> Xhat=Xh)
**> }
**>
**> ============================================
**>
**> Now lets take your original data and the first Y-value
**> in your list (namely Y = 1.251), and suppose you want
**> a 95% confidence interval for X. The X-value corresponding
**> to Y which you would get by regressing x (conc) on y (abs)
**> is X = 131.3813 so use this as a "starting value".
**>
**> So now run this function with x<-conc, y<-abs, and these values
**> of X and Y:
**>
**> R.calib(x,y,131.3813,1.251)
**>
**> You get a long list of stuff, amongst which
**>
**> $PF
**> [1] 0.02711878
**>
**> and
**>
**> $Xhat
**> [1] 131.2771
**>
**> So now you know that Xhat (the MLE of X for that Y) = 131.2771
**> and the F-ratio probability is 0.027...
**>
*

*****> You want to push $PF upwards till it reaches 0.05, so work
*****> *outwards* in the X-value:

WRONG!! Till it reaches ***0.95***

R.calib(x,y,125.0000,1.251)$PF

[1] 0.9301972

...

R.calib(x,y,124.3510,1.251)$PF

[1] 0.949989

> and you're there in that direction. Now go back to the MLE

*> and work out in the other direction:
**>
**> R.calib(x,y,131.2771,1.251)$PF
**> [1] 1.987305e-06
*

...

R.calib(x,y,138.0647,1.251)$PF

[1] 0.95

> and again you're there.

*>
**> So now you have the MLE Xhat = 131.2771, and the two
*

****> limits of a 95% confidence interval (131.0847, 131.4693)
**WRONG!!!
**

limits of a confidence interval (124.3510, 138.0647)

> for X, corresponding to the given value 1.251 of Y.

As it happens, these seem to correspond very closely to what you would get by reading "predict" in reverse:

*> plot(x,y)
*

> plm<-lm(y~x)

> min(x)

[1] 100

> max(x)

[1] 280

*> points((131.2771),(1.251),col="red",pch="+") #The MLE of X
**> lines(c(124.3506,138.0647),c(1.251,1.251),col="red") # The above CI
**> newx<-data.frame(x=(100:280))
**> predlm<-predict(plm,newdata=newx,interval="prediction")
**> lines(newx$x,predlm[,"fit"],col="green")
*

> lines(newx$x,predlm[,"lwr"],col="blue")

> lines(newx$x,predlm[,"upr"],col="blue")

which is what I thought would happen in the first place, given the quality of your data.

Sorry for any inconvenience due to the above error. Ted.

E-Mail: (Ted Harding) <Ted.Harding@nessie.mcc.ac.uk> Fax-to-email: +44 (0)870 094 0861

Date: 20-Apr-05 Time: 08:58:37 ------------------------------ XFMail ------------------------------ ______________________________________________R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Wed Apr 20 18:58:32 2005

*
This archive was generated by hypermail 2.1.8
: Fri 03 Mar 2006 - 03:31:17 EST
*