9

I would like to get the standard error on a prediction. Using R glm, I can get the SE of the fit for a specific prediction:

mod <- glm(y~wa_WSI, data=mydata, family=gaussian(link="identity"))
predict.glm(mod,newdata=newdata, type="response", se.fit=T)

But when I compare the predictions with the actual values, this number seems way too small. I found a formula for "standard error of the estimate" which is $\sqrt{s/(n-p)}$ where $s$ is the sum of the squared residuals, $n$ is the number of data points, and $p$ is the number of terms in the regression. This gives me a much larger result, but is not for a single prediction.

My question is, is the SE formula above the formula I should use and is there some way to get it from the value R gives me for se.fit so that it is specific for a particular prediction?

gung - Reinstate Monica
  • 132,789
  • 81
  • 357
  • 650
John
  • 425
  • 1
  • 4
  • 8
  • @GregSnow has provided a helpful answer. For more information about prediction intervals, it may help you to read this: [linear-regression-prediction-interval](http://stats.stackexchange.com/questions/33433//33642#33642). – gung - Reinstate Monica Feb 14 '14 at 20:38
  • 2
    The se.fit that predict.glm produces is a standard error for the *mean* prediction. For some GLMs it's meaningful to talk about a prediction interval (e.g. for the normal and Gamma), and a standard error for a future observation, but even in the cases where it makes sense, the problem - while easy for the normal - is difficult in the general case. You can do (for example) an asymptotic simulation -- simulate from $\hat\eta-\eta$ and then from $(Y|\eta-\hat\eta)$. – Glen_b Feb 14 '14 at 23:04
  • If you do want to compute the standard error on your predictions using `se.fit`, you should be able to do so as follows: `sqrt(predict(mod, newdata, se.fit = TRUE)$se.fit^2 + predict(mod, newdata, se.fit = TRUE)$residual.scale^2)`. – DataProphets Apr 19 '20 at 16:06

2 Answers2

7

I would like to throw in a comment for non-normal distributions and non-identity link functions. se.fit=T yields standard errors of the prediction, i.e. a measure of uncertainty for the predicted value. This prediction, by one of the Central Value Theorems, can be assumed to be normally distributed at the link scale, and hence its standard error can be given as the standard deviation of a normal distribution.

When using type="response", the prediction is back-transformed with the anti-link function (e.g. plogis for the logit-link). Using type="response"and se.fit=T yields non-sensical values, as it only returns one set of standard errors at the response scale. As the link-function is non-linear, the symmetric errors at the link scale must be asymmetric at the response scale. Thus, we can choose type="response" or se.fit=T, but not both when using non-identity link functions. (I don't understand why predict.glm has not been programmed to throw an error in this case.)

Julien
  • 3
  • 4
Carsten
  • 341
  • 2
  • 3
4

It is hard to answer without knowing more about what mod is. That is why we suggest a reproducible example.

If mod is a glm fit with a 'gaussian' family (the default) then it is just a linear model and you can use predict.lm instead which has the interval argument that can be set to "prediction" to compute prediction intervals.

If mod is a glm fit with a non-Gaussian family then the concept of a standard error of prediction may not even make sense (what is the prediction interval when the predictions are all TRUE/FALSE?).

If you can give more detail (a reproducible example and a clear statement of what you want) then we will have a better chance of giving a useful answer.

Greg Snow
  • 46,563
  • 2
  • 90
  • 159
  • 1
    I am so happy you pointed this simple thing out to me!! Yes my model was actually linear so I changed it to use lm and I have the CIs. I am still a little confused about what the se.fit actually is.. – John Feb 14 '14 at 19:59