多项式回归的置信区间

tra*_*syx 5 regression r confidence-interval polynomials

我对R和统计有一点疑问。

我使用最大似然法对模型进行拟合,该模型为我提供了以下系数以及各自的标准误差(除了其他参数估计值):

    ParamIndex   Estimate     SE        
1         a0  0.2135187 0.02990105  
2         a1  1.1343072 0.26123775  
3         a2 -1.0000000 0.25552696  
Run Code Online (Sandbox Code Playgroud)

从我可以画出的曲线:

 y= 0.2135187 + 1.1343072 * x - 1 * I(x^2)
Run Code Online (Sandbox Code Playgroud)

但是从那以后,我现在必须计算该曲线周围的置信区间,而且我不清楚如何做到这一点。

显然,我应该使用传播或错误/不确定性,但是我发现的方法需要原始数据,或者不仅仅是多项式公式。

当估计的SE已知为R时,有什么方法可以计算曲线的CI?

谢谢您的帮助。


编辑:

所以,现在,我有了函数的协方差表(v)vcov

                 a0           a1           a2
    a0  0.000894073 -0.003622614  0.002874075
    a1 -0.003622614  0.068245163 -0.065114661
    a2  0.002874075 -0.065114661  0.065294027
Run Code Online (Sandbox Code Playgroud)

n = 279

李哲源*_*李哲源 6

您现在没有足够的信息。要计算拟合曲线的置信区间,需要三个系数的完整方差-协方差矩阵,但现在您只有该矩阵的对角线条目。

如果您已拟合正交多项式,则方差-协方差矩阵是对角矩阵,具有相同的对角元素。这当然不是你的情况,因为:

  • 您显示的标准错误彼此不同;
  • 您已明确使用原始多项式表示法:x + I(x ^ 2)

但我发现的方法需要原始数据

它不是用于拟合模型的“原始数据”。这是您想要生成置信带的“新数据”。但是,您确实需要知道用于拟合模型的数据数量,例如n,因为这是导出剩余自由度所必需的。在有 3 个系数的情况下,这个自由度是n - 3

一旦你拥有:

  • 完整的方差-协方差矩阵,比方说V
  • n,用于模型拟合的数据数量;
  • x给出在何处产生置信带的点向量,

您可以首先从以下位置获取预测标准误差:

X <- cbind(1, x, x ^ 2)    ## prediction matrix
e <- sqrt( rowSums(X * (X %*% V)) )    ## prediction standard error
Run Code Online (Sandbox Code Playgroud)

您知道如何从拟合多项式公式中获得预测平均值,对吧?假设平均值为mu,现在对于 95%-CI,使用

## residual degree of freedom: n - 3
mu + e * qt(0.025, n - 3)  ## lower bound
mu - e * qt(0.025, n - 3)  ## upper bound
Run Code Online (Sandbox Code Playgroud)

完整的理论位于Predict.lm() 如何计算置信区间和预测区间?


更新

根据您提供的协方差矩阵,现在可以产生一些结果和数字。

V <- structure(c(0.000894073, -0.003622614, 0.002874075, -0.003622614, 
0.068245163, -0.065114661, 0.002874075, -0.065114661, 0.065294027
), .Dim = c(3L, 3L), .Dimnames = list(c("a0", "a1", "a2"), c("a0", 
"a1", "a2")))
Run Code Online (Sandbox Code Playgroud)

假设我们想在 处生成 CI x = seq(-5, 5, by = 0.2)

beta <- c(0.2135187, 1.1343072, -1.0000000)
x <- seq(-5, 5, by = 0.2)
X <- cbind(1, x, x ^ 2)
mu <- X %*% beta
e <- sqrt( rowSums(X * (X %*% V)) )
n <- 279
lo <- mu + e * qt(0.025, n - 3)
up <- mu - e * qt(0.025, n - 3)
matplot(x, cbind(mu, lo, up), type = "l", col = 1, lty = c(1,2,2))
Run Code Online (Sandbox Code Playgroud)

在此输入图像描述