R help - Nov 2023 - understanding predict.lm

Hello, All:


	  I am unable to manually replicate predict.lm, specifically comparing 
se.fit with (fit[,3]-fit[,2]): I think their ratio should be 
2*qnorm((1-level)/2), and that's not what I'm getting.


	  Consider the following slight modification of the first example in 
help('predict.lm'):


set.seed(1)
x <- rnorm(15)
y <- x + rnorm(15)
predict(lm(y ~ x))
new <- data.frame(x = seq(-3, 3, 0.5))
predict(lm(y ~ x), new, se.fit = TRUE)
pred.w.plim <- predict(lm(y ~ x), new, interval = "prediction",
                        se.fit = TRUE)
pred.w.clim <- predict(lm(y ~ x), new, interval = "confidence",
                        se.fit = TRUE)

(z.confInt <- with(pred.w.clim, (fit[,3]-fit[,2])/se.fit))
pnorm(-z.confInt/2)

s.pred <- sqrt(with(pred.w.plim,
                     se.fit^2+residual.scale^2))
(z.predInt <- with(pred.w.plim, (fit[,3]-fit[,2])/s.pred))
pnorm(-z.predInt/2)


	  ** This gives me 0.01537207. I do not understand why it's not 0.025 
with level = 0.95.


	  Can someone help me understand this?
	  Thanks,
	  Spencer Graves

Dear Spencer,

You need the t distribution with correct df, not the standard-normal 
distribution:

 > pt(-z.confInt/2, df=13)
     1     2     3     4     5     6     7     8     9    10    11
0.025 0.025 0.025 0.025 0.025 0.025 0.025 0.025 0.025 0.025 0.025
    12    13
0.025 0.025

 > pt(-z.predInt/2, df=13)
     1     2     3     4     5     6     7     8     9    10    11
0.025 0.025 0.025 0.025 0.025 0.025 0.025 0.025 0.025 0.025 0.025
    12    13
0.025 0.025

I hope this helps,
  John
-- 
John Fox, Professor Emeritus
McMaster University
Hamilton, Ontario, Canada
web: https://www.john-fox.ca/

On 2023-11-06 4:16 p.m., Spencer Graves wrote:
> Caution: External email.
> 
> 
> Hello, All:
> 
> 
>  ???????? I am unable to manually replicate predict.lm, specifically 
> comparing
> se.fit with (fit[,3]-fit[,2]): I think their ratio should be
> 2*qnorm((1-level)/2), and that's not what I'm getting.
> 
> 
>  ???????? Consider the following slight modification of the first 
> example in
> help('predict.lm'):
> 
> 
> set.seed(1)
> x <- rnorm(15)
> y <- x + rnorm(15)
> predict(lm(y ~ x))
> new <- data.frame(x = seq(-3, 3, 0.5))
> predict(lm(y ~ x), new, se.fit = TRUE)
> pred.w.plim <- predict(lm(y ~ x), new, interval =
"prediction",
>  ?????????????????????? se.fit = TRUE)
> pred.w.clim <- predict(lm(y ~ x), new, interval =
"confidence",
>  ?????????????????????? se.fit = TRUE)
> 
> (z.confInt <- with(pred.w.clim, (fit[,3]-fit[,2])/se.fit))
> pnorm(-z.confInt/2)
> 
> s.pred <- sqrt(with(pred.w.plim,
>  ??????????????????? se.fit^2+residual.scale^2))
> (z.predInt <- with(pred.w.plim, (fit[,3]-fit[,2])/s.pred))
> pnorm(-z.predInt/2)
> 
> 
>  ???????? ** This gives me 0.01537207. I do not understand why it's not
> 0.025
> with level = 0.95.
> 
> 
>  ???????? Can someone help me understand this?
>  ???????? Thanks,
>  ???????? Spencer Graves
> 
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide 
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

Doh! Thanks very much. sg


On 11/6/23 5:17 PM, John Fox wrote:
> Dear Spencer,
> 
> You need the t distribution with correct df, not the standard-normal 
> distribution:
> 
>  > pt(-z.confInt/2, df=13)
>  ??? 1???? 2???? 3???? 4???? 5???? 6???? 7???? 8???? 9??? 10??? 11
> 0.025 0.025 0.025 0.025 0.025 0.025 0.025 0.025 0.025 0.025 0.025
>  ?? 12??? 13
> 0.025 0.025
> 
>  > pt(-z.predInt/2, df=13)
>  ??? 1???? 2???? 3???? 4???? 5???? 6???? 7???? 8???? 9??? 10??? 11
> 0.025 0.025 0.025 0.025 0.025 0.025 0.025 0.025 0.025 0.025 0.025
>  ?? 12??? 13
> 0.025 0.025
> 
> I hope this helps,
>  ?John