Re: [R] Help with predict.lm

From: Mike White <mikewhite.diu_at_tiscali.co.uk>
Date: Wed 20 Apr 2005 - 20:16:12 EST

Ted
Thank you for giving me so much of your time to provide an explanation of the likelihood ratio approach to the calibration problem. It has certainly provided me with a new insight into this problem. I will check out the references you mentioned to get a better understanding of the statistics. Out of interest I have also tried the function provided by David Lucy back in March 2002 with a correction to the code and addition of the t-value. The results are very close to those obtained by the likelihood ratio method, although the output needs tidying up.

####################

# R function to do classical calibration # val is a data.frame of the new indep variables to predict dep

calib <- function(indep, dep, val, prob=0.95) {

 # generate the model y on x and get out predicted values for x

        reg <- lm(dep~indep)
        xpre <- (val - coef(reg)[1])/coef(reg)[2]

        # generate a confidence
        yyat <- ((dep - predict(reg))^2)
        sigyyat <- ((sum(yyat)/(length(dep)-2))^0.5)
        term1 <- sigyyat/coef(reg)[2]
        sigxxbar <- sum((indep - mean(indep))^2)
        denom <- sigxxbar * ((coef(reg)[2])^2)
 t<-qt(p=(prob+(1-prob)/2), df=(length(dep)-2))
        conf <- t*abs((((1+1/length(dep))+(((val -
  mean(dep))^2)/denom))^0.5)*term1)

results<-list(Predicted=xpre, Conf.interval=c(xpre-conf, xpre+conf)) return(results)
}

conc<-seq(100, 280, 20) # mg/l
abs<-c(1.064, 1.177, 1.303, 1.414, 1.534, 1.642, 1.744, 1.852, 1.936, 2.046) # absorbance units

new.abs<-data.frame(abs=c(1.251, 1.324, 1.452))

calib(conc, abs, new.abs, prob=0.95)

$Predicted

       abs

1 131.2771
2 144.6649
3 168.1394

$Conf.interval
$Conf.interval$abs

[1] 124.4244 137.9334 161.5477

$Conf.interval$abs

[1] 138.1298 151.3964 174.7310

Thanks
Mike

> 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 23:25:33 2005

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