代码之家  ›  专栏  ›  技术社区  ›  IceCreamToucan

R中lm对象的预测函数

  •  1
  • IceCreamToucan  · 技术社区  · 8 年前

    为什么是 prediction_me prediction_R 不相等?我试图遵循引理5给出的公式 here . 是否 predict 函数使用不同的公式,我的计算是否有错误,或者只是舍入错误?(两者非常接近)

    set.seed(100)
    # genrate data
    x    <- rnorm(100, 10)
    y    <- 3 + x + rnorm(100, 5)
    data <- data.frame(x = x, y = y)
    # fit model
    mod  <- lm(y ~ x, data = data)
    
    # new observation
    data2 <- data.frame(x = rnorm(5, 10))
    
    # prediction for new observation
    d    <- as.matrix(cbind(1, data[,-2]))
    d2   <- as.matrix(cbind(1, data2))
    fit  <- d2 %*% mod$coefficients 
    t    <- qt(1 - .025, mod$df.residual)
    s    <- summary(mod)$sigma
    half <- as.vector(t*s*sqrt(1 + d2%*%solve(t(d)%*%d, t(d2))))
    
    prediction_me <- cbind(fit, fit - half, fit + half)
    
    prediction_R <- predict(mod, newdata = data2, interval = 'prediction')
    
    
    prediction_me
    prediction_R
    
    1 回复  |  直到 8 年前
        1
  •  2
  •   Julius Vainora    8 年前

    您当前的代码几乎没有问题。只需注意引理5中的公式是针对一个新观察到的 x . 因此, half 不仅包含相关方差,还包含协方差,而您只需要前一个方差。因此 as.vector 应替换为 diag :

    half <- diag(t * s * sqrt(1 + d2 %*% solve(t(d) %*%d , t(d2))))
    prediction_me <- cbind(fit, fit - half, fit + half)
    prediction_R <- predict(mod, newdata = data2, interval = 'prediction')
    
    range(prediction_me - prediction_R)
    # [1] 0 0