[Rd] Problems with predict.lm: incorrect SE estimate (PR#7772)

[Peter Dalgaard]

murdoch at stats.uwo.ca writes:

> On Mon,  4 Apr 2005 18:01:05 +0200 (CEST), msa at biostat.mgh.harvard.edu
> wrote :
> 
> >Full_Name: Marek Ancukiewicz
> >Version: 2.01
> >OS: Linux
> >Submission from: (NULL) (132.183.12.87)
> >
> >
> >It seems that the the standard error of prediction of the linear regression,
> >caclulated with predict.lm is incorrect. Consider the following example where
> >the standard error is first calculated with predict.lm, then using delta
> >method. and finally, using the formula rms*sqrt(1+1/n+(xp-x0)^2/Sxx). 
> 
> Your formula is incorrect.  You've got the formula for the so called
> "prediction error" (i.e. the stddev of the difference between the
> prediction and a new observation) rather than the "standard error"
> (i.e. the stddev of the prediction).

And:

>  print(predict(a,new=data.frame(x=xp),interval="pred"))
          fit      lwr      upr
[1,] 3.009091 2.794523 3.223659
>  3.009091  + qt(.975,8)*0.09304758
[1] 3.223659
>  3.009091  - qt(.975,8)*0.09304758
[1] 2.794523

so reading the help page might have given a clue that the authors knew
what they were doing....

The help page text could be improved, though. Will do.

> >$fit
> >[1] 3.009091
> >
> >$se.fit
> >[1] 0.0359752
> >
> >$df
> >[1] 8
> >
> >$residual.scale
> >[1] 0.08581163
> >
> >> print(se.delta.method <- sqrt(s[1,1]+xp^2*s[2,2]+2*xp*s[1,2] + rms^2))
> >[1] 0.09304758
> >> print(se.ss.formula <- rms*sqrt(1+1/n+(xp-xm)^2/sum((x-xm)^2)))
> >[1] 0.09304758

-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)             FAX: (+45) 35327907