首页 文章

多项式回归的置信区间

提问于
浏览
4

我和R和统计数据有点问题 .

我使用最大似然法拟合了一个模型,他给出了以下系数及其各自的标准误差(以及其他参数估计值):

ParamIndex   Estimate     SE        
1         a0  0.2135187 0.02990105  
2         a1  1.1343072 0.26123775  
3         a2 -1.0000000 0.25552696

从我可以画出我的曲线:

y= 0.2135187 + 1.1343072 * x - 1 * I(x^2)

但是从那时起,我现在要计算这条曲线周围的置信区间,我不清楚如何做到这一点 .

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

当用R知道估计值的SE时,有没有任何方法来计算我的曲线的CI?

谢谢您的帮助 .


Edit:

所以,现在,我使用函数 vcov 获得协方差表(v):

a0           a1           a2
    a0  0.000894073 -0.003622614  0.002874075
    a1 -0.003622614  0.068245163 -0.065114661
    a2  0.002874075 -0.065114661  0.065294027

n = 279 .

1 回答

  • 3

    You don't have enough information right now. 要计算拟合曲线的置信区间,需要使用 a complete variance-covariance matrix 作为三个系数,但现在只有该矩阵的对角线条目 .

    如果已经拟合了正交多项式,则方差 - 协方差矩阵是对角线的,具有相同的对角线元素 . 这肯定不是你的情况,因为:

    您显示的

    • 标准错误彼此不同;

    • 您已明确使用原始多项式表示法: x + I(x ^ 2)

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

    它不是用于拟合模型的"raw data" . 它是"new data"你想要产生信心带 . 但是,您确实需要知道用于拟合模型的数据的数量,例如 n ,因为这是获得剩余自由度所必需的 . 在你的3个系数的情况下,这个自由度是 n - 3 .

    一旦你有:

    • 完全方差 - 协方差矩阵,让我们说 V ;

    • n ,用于模型拟合的数据数量;

    • x 的向量给出了产生置信带的位置,

    您可以先从以下方面获得预测标准误差:

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

    你知道如何从拟合的多项式公式得到预测均值吗?假设平均值为 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
    

    一个完整的理论是How does predict.lm() compute confidence interval and prediction interval?


    Update

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

    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")))
    

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

    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))
    

    enter image description here

相关问题