我跑回了一个回归:
CopierDataRegression <- lm(V1~V2, data=CopierData1)
我的任务是获得一个
- 90% confidence interval 为
V2=6
和V2=6
的平均响应
当V2=6
时 - 90% prediction interval .
我使用了以下代码:
X6 <- data.frame(V2=6)
predict(CopierDataRegression, X6, se.fit=TRUE, interval="confidence", level=0.90)
predict(CopierDataRegression, X6, se.fit=TRUE, interval="prediction", level=0.90)
我得到 (87.3, 91.9)
和 (74.5, 104.8)
似乎是正确的,因为PI应该更宽 .
两者的输出也包括 se.fit = 1.39
,这是相同的 . I don't understand what this standard error is. Shouldn't the standard error be larger for the PI vs. the CI? How do I find these two different standard errors in R?
数据:
CopierData1 <- structure(list(V1 = c(20L, 60L, 46L, 41L, 12L, 137L, 68L, 89L,
4L, 32L, 144L, 156L, 93L, 36L, 72L, 100L, 105L, 131L, 127L, 57L,
66L, 101L, 109L, 74L, 134L, 112L, 18L, 73L, 111L, 96L, 123L,
90L, 20L, 28L, 3L, 57L, 86L, 132L, 112L, 27L, 131L, 34L, 27L,
61L, 77L), V2 = c(2L, 4L, 3L, 2L, 1L, 10L, 5L, 5L, 1L, 2L, 9L,
10L, 6L, 3L, 4L, 8L, 7L, 8L, 10L, 4L, 5L, 7L, 7L, 5L, 9L, 7L,
2L, 5L, 7L, 6L, 8L, 5L, 2L, 2L, 1L, 4L, 5L, 9L, 7L, 1L, 9L, 2L,
2L, 4L, 5L)), .Names = c("V1", "V2"),
class = "data.frame", row.names = c(NA, -45L))
2 回答
指定
interval
和level
参数时,predict.lm
可以返回置信区间(CI)或预测区间(PI) . 此答案显示如何在不设置这些参数的情况下获取CI和PI . 有两种方法:使用
predict.lm
的中期结果;从头开始做一切 .
了解如何使用这两种方式可以让您全面了解预测过程 .
请注意,我们只会涵盖
predict.lm
的type = "response"
(默认)大小写 . 讨论type = "terms"
超出了这个答案的范围 .设置
我在这里收集你的代码,以帮助其他读者复制,粘贴和运行 . 我还更改变量名称,以便它们具有更清晰的含义 . 另外,我扩展
newdat
以包含多行,以显示我们的计算是"vectorized" .以下是
predict.lm
的输出,稍后将与我们的手动计算进行比较 .使用来自predict.lm的中间阶段结果
z$se.fit
是预测平均值z$fit
的标准误差,用于构造z$fit
的CI . 我们还需要具有自由度的分布的分位数z$df
.我们看到这与
predict.lm(, interval = "confidence")
一致 .PI比CI更宽,因为它考虑了剩余方差:
请注意,这是逐点定义的 . 对于非加权线性回归(如在您的示例中),残差方差在任何地方都相等(称为同方差性),并且它是
z$residual.scale ^ 2
. 因此,PI的标准误差是PI构建为
我们看到这与
predict.lm(, interval = "prediction")
一致 .remark
如果你有一个权重线性回归,那么事情会更复杂,其中残差方差在任何地方都不相等,因此应该加权.756734_ . 为拟合值构造PI更容易(也就是说,在
predict.lm
中使用type = "prediction"
时不设置newdata
),因为权重是已知的(使用lm
时必须通过weight
参数提供) . 对于样本外预测(即,您将newdata
传递给predict.lm
),predict.lm
希望您告诉它应如何对残差方差进行加权 . 您需要在predict.lm
中使用参数pred.var
或weights
,否则您会收到来自predict.lm
的警告,抱怨构建PI的信息不足 . 以下引用自?predict.lm
:请注意,CI的构造不受回归类型的影响 .
从头开始做一切
基本上我们想知道如何在
z
中获得fit
,se.fit
,df
和residual.scale
.预测平均值可以通过矩阵向量乘法
Xp %*% b
来计算,其中Xp
是线性预测矩阵,b
是回归系数向量 .我们看到这与
z$fit
一致 .yh
的方差 - 协方差是Xp %*% V %*% t(Xp)
,其中V
是b
的方差 - 协方差矩阵,可以通过不需要
yh
的完全方差 - 协方差矩阵来计算逐点CI或PI . 我们只需要它的主对角线 . 因此,我们可以通过提高效率来实现diag(Xp %*% V %*% t(Xp))
在拟合模型中可以轻松获得剩余自由度:
最后,要计算残差方差,请使用Pearson估算器:
remark
请注意,在加权回归的情况下,
sig2
应计算为附录:一个模仿predict.lm的自编函数
"Do everything from scratch"中的代码已完整地组织到此问答中的函数
lm_predict
:linear model with lm: how to get prediction variance of sum of predicted values .我不知道是否有一种快速的方法来提取预测间隔的标准误差,但是你可以随时回溯SE的间隔(即使它不是超级优雅的方法):
请注意,CI SE与
se.fit
中的值相同 .