単回帰分析のRのコードをガチで実装してみる

single_regression.R

single_regr <- function(mdl, xh, p){
x <- mdl$model[,2]; y <- mdl$model[,1]; e <- mdl$resid; n<- length(x)
a <- mdl$coef[1]; b <- mdl$coef[2]
yh <- a + b*xh
sgm <- sqrt(sum(e*e)/(n-2)); sx <- var(x)*(n-1)
ms <- sqrt(1+(1/n) + (xh-mean(x))^2/sx)*sgm
k <- (1-p)/2; tk <- qt((1-k),df=(n-2))
yL <- yh - tk*ms; yU <- yh + tk*ms
xmin <- min(xh,x); xmax <- max(xh,x)
ymin <- min(yL,y); ymax <- max(yU,y)
title <- paste("single-regression Confidence interval: ",p)
plot(x, y, xlim=c(xmin, xmax), ylim=c(ymin, ymax), main=title)
lines(xh,yh,lty=1); lines(xh,yL, lty=2, col="red")
lines(xh, yU, lty=2, col="blue")
return(list(xh=xh, yL=yL, yh=yh, yU=yU))
}