Я пытаюсь вычислить доверительный интервал вокруг линейной регрессии в R. Я знаю, что predict может сделать это в большинстве случаев, но я хочу решение на основе уравнения; отчасти потому, что не все регрессионные модели (например, из пакета deming) совместимы с predict, а отчасти потому, что я хочу это понять (отсюда также нет решения ggplot или подобного).
I зашли довольно далеко, используя уравнения, объясненные здесь и здесь .
Что я достиг: Я успешно вычислил t -значение MSE и стандартная ошибка регрессии. Я очень уверен, что все уравнения в моей реализации в значительной степени корректны, потому что, когда я использую их для диапазона значений x , на котором была рассчитана регрессионная модель, они идеально соответствуют тому, что predict возвращает (см. левый график в коде).
Где это усложняется: Проблемы начинаются, когда я пытаюсь экстраполировать; или другими словами, когда я хочу вычислить доверительный интервал за пределами диапазона x -значений, для которых у меня есть данные. Рассчитанный доверительный интервал все еще корректен , но он смещен . Точно, как вы можете видеть на правом графике ниже, вы должны переместить доверительный интервал вдоль оси x на разницу средних значений исходных данных и данных, использованных для экстраполяции. Аналогично, на оси y полоса достоверности должна быть смещена разницей модели для соответствующих средних значений x . Вы видите вычисления под #Define displacement vector в коде ниже, если это трудно понять.
Я смутно понимаю, почему происходит это смещение, учитывая уравнения для se. Но мне интересно, есть ли лучшее решение проблемы, чем замена полосы доверия, которую я реализовал сейчас (также потому, что из-за необходимости смещения мой код теперь не действительно рассчитывает полосу доверия через предполагаемый интервал). Я был бы очень признателен, если бы кто-нибудь мог помочь мне с этим кодом.
#Create data
Dat<-as.data.frame(matrix(c(1, 1, 1, 4, 4, 4, 7, 7, 7, 10, 10, 10, 2.1, 2.3, 2.2, 3.5, 3.1, 3.2, 4.2, 5.0, 4.8, 6.1, 6.6, 6.2), 12, 2))
colnames(Dat)<-c("X", "Y")
#Create linear model
mod<-lm(Y ~ X, data=Dat)
#Use predict to calculate confidence band for comparison
Pred<-predict(mod, newdata=data.frame(X=0:30), interval="confidence")
#Calculate confidence band according to equations
#https://stattrek.com/regression/slope-confidence-interval.aspx
#https://library2.lincoln.ac.nz/documents/Analysing-the-Variance.pdf
##Gather constants
n<-nrow(Dat)
##Define prediction values
Pred.vals<-list()
Pred.vals$S1<-seq(from=1, to=10, by=0.5)
Pred.vals$S2<-seq(from=0, to=30, by=0.5)
Pred.vals$S1.fitted<-coef(mod)[2]*Pred.vals$S1+coef(mod)[1]
Pred.vals$S2.fitted<-coef(mod)[2]*Pred.vals$S2+coef(mod)[1]
##Calculate t-value
t.val<-qt(p=1-((1-0.95)/2), df=n-2)
##Calculate MSE
mse<-sqrt(sum((Dat[,"Y"]-mod$fitted.values)^2)/(n-2))
##Calculate standard error of fit: two versions, both work, but se2 is displaced
se1<-mse*sqrt((1/n)+(Pred.vals$S1-mean(Pred.vals$S1))^2/sum((Dat[,"X"]-mean(Dat[,"X"]))^2))
se2<-mse*sqrt((1/n)+(Pred.vals$S2-mean(Pred.vals$S2))^2/sum((Dat[,"X"]-mean(Dat[,"X"]))^2))
#Define displacement vector
X.Mean<-list()
X.Mean$Original<-mean(Dat[,"X"])
X.Mean$New<-mean(Pred.vals$S2)
X.Mean$X.Diff<-X.Mean$Original-X.Mean$New
X.Mean$Y.Diff<-(coef(mod)[2]*X.Mean$Original+coef(mod)[1])-(coef(mod)[2]*X.Mean$New+coef(mod)[1])
#Calculate confidence band
slope.upper1<-Pred.vals$S1.fitted+t.val*se1
slope.lower1<-Pred.vals$S1.fitted-t.val*se1
slope.upper2<-Pred.vals$S2.fitted+t.val*se2
slope.lower2<-Pred.vals$S2.fitted-t.val*se2
#Plot and compare
win.graph(20, 10, 10)
layout(matrix(1:2, 1, 2))
##Small plot
plot(Dat[,"X"], Dat[,"Y"], xlim=c(0, 11), ylim=c(2, 7))
curve(coef(mod)[2]*x+coef(mod)[1], col="grey50", lwd=1, add=TRUE)
##Confidence interval from "predict"
lines(0:30, Pred[,"lwr"], col="cornflowerblue", lty=2)
lines(0:30, Pred[,"upr"], col="cornflowerblue", lty=2)
##Confidence intervals from equations
lines(Pred.vals$S1, slope.upper1, col="darkgreen", lwd=2, lty=2)
lines(Pred.vals$S1, slope.lower1, col="darkgreen", lwd=2, lty=2)
legend("topleft", col=c("grey50", "cornflowerblue", "darkgreen"), lwd=c(1, 1, 2), lty=c(1, 2, 2), legend=c("Regression line", "Confidence from 'predict'", "Confidence from equations"))
##Large plot
plot(Dat[,"X"], Dat[,"Y"], xlim=c(0, 30), ylim=c(2, 15))
curve(coef(mod)[2]*x+coef(mod)[1], col="grey50", lwd=1, add=TRUE)
##Confidence interval from "predict"
lines(0:30, Pred[,"lwr"], col="cornflowerblue", lty=2)
lines(0:30, Pred[,"upr"], col="cornflowerblue", lty=2)
##Confidence intervals from equations
#lines(Pred.vals$S1, slope.upper1, col="darkgreen", lwd=2, lty=2)
#lines(Pred.vals$S1, slope.lower1, col="darkgreen", lwd=2, lty=2)
lines(Pred.vals$S2, slope.upper2, col="firebrick", lty=3)
lines(Pred.vals$S2, slope.lower2, col="firebrick", lty=3)
lines(Pred.vals$S2+X.Mean$X.Diff, slope.upper2+X.Mean$Y.Diff, col="darkgreen", lwd=2, lty=3)
lines(Pred.vals$S2+X.Mean$X.Diff, slope.lower2+X.Mean$Y.Diff, col="darkgreen", lwd=2, lty=3)
legend("topleft", col=c("grey50", "cornflowerblue", "firebrick", "darkgreen"), lwd=c(1, 1, 1, 2), lty=c(1, 2, 3, 3), legend=c("Regression line", "Confidence from 'predict'", "Confidence from equations", "Confidence from equations (displaced)"))