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) <[EMAIL PROTECTED]>
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
