r 函数:多元线性回归预测估计和区间(用户定义的函数)

r function: multiple linear regression prediction estimate and interval (user-defined function)

我正在使用 r 中的用户定义函数来计算预测估计值和 95% 线性回归的区间。我有一个函数可以复制 predict.lm() 函数拟合和间隔。但是,当应用于多元线性回归时,小数点后第三位略有不同,我无法解释原因。我没有扎实的理论数学背景,所以我使用了这个网站和解释的公式来整合到我的函数中:https://daviddalpiaz.github.io/appliedstats/multiple-linear-regression.html

我的函数是否存在错误,或者由于舍入误差或其他边际误差而造成的细微差异?下面是功能代码以及我如何应用它来测试它:

predict.reg.95 <- function(lm.model,newdata) {
  if (!inherits(lm.model, "lm")){warning("object is not a lm() model")}
  else{
    n<-length(lm.model$residuals)
    beta <- lm.model$coefficients
    sy<-sigma(lm.model)
    s2x<-var(lm.model$model[,2])
    t.alpha.demi<- qt(0.975, df=n-2)
      if (length(beta)-1==length(newdata)) {
        y.pred <- beta[1]+sum(beta[-1]*newdata)
        x0<-c(1,newdata)
        X<-cbind(c(rep(1,n)),lm.model$model[,-1])
        y.pred.interval.upp<-y.pred+t.alpha.demi*sy*sqrt(1+x0%*%solve(t(X)%*%as.matrix(X))%*%x0)
        y.pred.interval.low<-y.pred-t.alpha.demi*sy*sqrt(1+x0%*%solve(t(X)%*%as.matrix(X))%*%x0)
        fit<-c(y.pred)
        upr<-c(y.pred.interval.upp)
        lwr<-c(y.pred.interval.low)
        output<-cbind(lwr,fit,upr)
        print("Below, you will find the predicted estimate (fit) with the given values of the explanatory variables and the associated prediction interval (lwr,upr)")
        print(output)
      } 
      else {
      print("the length of the chosen explanatory variables vector isn't the same length as the number of explanatory variables of your lm() model")
     }
    }
  }
#1st test
df<-data.frame(x1=c(sample.int(100,50, replace=T)),y=c(sample.int(200,50, replace=T)),x2=c(sample.int(20,50, replace=T)))
lm.model<-lm(y~x1+x2,data=df)
x1<-c(12)
x2<-c(32)
newdata<-as.data.frame(cbind(x1,x2))
new<-c(12,32)

predict.reg.95(lm.model,new)
predict(lm.model, newdata, level=0.95,interval="prediction") #slight difference at the third decimal for the prediction interval between functions

#2nd test
data(Seatbelts)
lm.model<-lm(DriversKilled~kms+drivers,data=Seatbelts)
kms<-c(10000)
drivers<-c(2000)
newdata<-as.data.frame(cbind(kms,drivers))
new<-c(10000,2000)

predict.reg.95(lm.model,new)
predict(lm.model, newdata, level=0.95,interval="prediction") #slight difference at the third decimal for the prediction interval between functions

希望有解决办法或者问题不大,功能可以原样使用

恭敬地,

西里尔 S

更新: 我发现了问题,它在 t 值 t.alpha.demi<- qt(0.975, df=n-2) 处,这解释了为什么它与单一线性回归没有区别但与多重线性回归有区别。

我改成了t.alpha.demi<- qt(0.975, df=n-length(beta))

这是我的一个错误。问候, 西里尔 S