Help with predict.lm

On 19-Apr-05 Mike White wrote:

Hi
I have measured the UV absorbance (abs) of 10 solutions
of a substance at known concentrations (conc) and have
used a linear model to plot a calibration graph with
confidence limits.  I now want to predict the concentration
of solutions with UV absorbance results given in the 
new.abs data.frame, however predict.lm only appears to work
for new "conc" variables not new "abs" variables.
[...]
    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
    lm.calibration<-lm(abs ~ conc)
    pred.w.plim <- predict(lm.calibration,  interval="prediction")
    pred.w.clim <- predict(lm.calibration,  interval="confidence")
    matplot(conc, cbind(pred.w.clim, pred.w.plim[,-1]),
       lty=c(1,2,2,3,3), type="l", ylab="abs", xlab= "conc mg/l")
    points(conc, abs, pch=21, col="blue")

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

[...]

The correct approach has in the past been the subject of
at times quite controversial discussion, under the title
indeed of "The Calibration Problem". Nowadays this problem
would be approached by making the concentrations to be
"predicted" additional unknown parameters, and evaluating
likelihood ratios for possible values of these.

I don't have time at the moment to go into this approach,
but will try to write something later.

Here is the basic theory, and an associated method, for
the case of estimating 1 value of X ("conc") corresponding
to 1 given value of Y ("abs"), when data x1,...,xn and
y1,...yn have been obtained (your data "conc" and "abs"
above, respectively).

The model

   y = a + b*x + e

is assumed, where the x-values are given (as typically in
a calibration experiment -- e.g. measuring standard test
samples of known y-value), and the y-values are measured
according to the above, with errors e assumed distributed
as N(0,s^2).

For simplicity, centre the x and y observations at their
means, so that Sum[xi]=0 and Sum[yi]=0.

Let an observation Y of y be given.
To infer the corresponding value of X.
(X is measured from the mean of the xi, Y from the mean of the yi).

Let X play the role of an additional parameter in the
problem. Then the likelihood function for (a,b,s,X) is

  (1/(s^(n+1))*exp(-(Sum[(y-a-b*x)^2] + (Y-a-b*X)^2)/s^2)

The MLEs of a#, b#, s#, X# of a, b, s, X are then given by

  a# = 0

  b# = (Sum[xy])/(Sum[x^2])

  X# = Y/b#

  (s#)^2 = (1/(n+1))Sum[(y - (b#)*x)^2

Now suppose that X is set at a fixed value (again denoted by X).
Then the MLEs a~, b~, s~ of a, b and s are now given by

  a~ = (Y*Sum[x^2] - X*Sum[x*y])/D

  b~ = ((n+1)*Sum[x*y] + n*X*Y)/D

  (s~)^2 = (Sum[(y - a~ - (b~)*x)^2] + (Y - a~ - (b~)*X)^2)/(n+1)

where

  D = (n+1)*Sum[x^2] + n*X^2

The likelihood ratio (profile likelihood for X) of the hypothesis
that X has a fixed value (X) versus the hypothesis that X might
be anything is therefore

  ((s#)/(s~))^(n+1)

which depends only on

  R(X) = (s#)^2/(s~)^2

where s~ (but not s#) depends on X.

Now

  (n-1)*((s~)^2 - (s#)^2)/(s#^2) = (n-2)*(1 - R(X))/R(X)

has the F distribution with 1 and (n-1) d.f., quite independent
of the true values of a, b and s^2, and large values of R(X)
correspond to small values of this "F ratio".

Hence the MLE of X is Y/B# and a confidence set for X at
a given confidence level P0 is the set of all X such that

  (n-2)*(1 - R(X))/R(X) < qf(1-P0,1,n-2)

(in the notation of R's F functions pf, qf, rf, etc.)

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:

  R.calib(x,y,131.4000,1.251)$PF
  [1] 0.03198323

  R.calib(x,y,131.4500,1.251)$PF
  [1] 0.0449835

  ...

  R.calib(x,y,131.4693,1.251)$PF
  [1] 0.04999851

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,131.2000,1.251)$PF
  [1] 0.02005908

  R.calib(x,y,131.1000,1.251)$PF
  [1] 0.04604698

  ...

  R.calib(x,y,131.0847,1.251)$PF
  [1] 0.0500181

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)
for X, corresponding to the given value 1.251 of Y.

I leave it to you (or any other interested parties)
to wrap this basic method in a proper routine which
will avoid the above groping (but at least it shows
explicitly what's going on).

As a refinement: In you original query, you put up
three values of Y: 1.251, 1.324, 1.452

Now you could, of course, simply apply the above
method to each one separately as above.

However, the same theoretical approach can be used
to obtain a joint confidence region for the three
corresponding X values jointly, and this is theoretically
more correct, since these are 3 extra parameters and
they will have a covariance structure. You would need
to work through the algebra for k such Y-values and the
corresponding X-values.

In practice this may not matter much -- probably not
at all for your data.

Note: The above analysis of the single-case situation
was exhibited in an address I gave to a modelling symposium
in 1985, subsequently published along with the other
presentations in a special number of The Statistician:

  Modelling: the classical approach
  The Statistician (1986) vol 35, pp. 115-134

NB that the equation for alpha-tilde (corresponding to a~ above)
in that article has a misprint ("Y" missing before Sum(x^2)).

This analysis had something in common with work by Phil Brown:

  P.J. Brown (1982) Multivariate Calibration
  JRSS Ser. B, vol 44, pp. 287-321

though Brown's approach did not start from the likelihood
ratio principle.

Hoping this helps!
Best wishes,
Ted.


--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at nessie.mcc.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 19-Apr-05                                       Time: 22:27:01
------------------------------ XFMail ------------------------------