2016-07-05 250 views
2

我有這樣的未來10年迴歸。產生良好的線性迴歸圖(擬合線,置信度/預測帶等)

date<-as.Date(c("2015-12-31", "2014-12-31", "2013-12-31", "2012-12-31")) 
value<-c(16348, 14136, 12733, 10737) 
#fit linear regression 
model<-lm(value~date) 
#build predict dataframe 
dfuture<-data.frame(date=seq(as.Date("2016-12-31"), by="1 year", length.out = 10)) 
#predict the futurne 
predict(model, dfuture, interval = "prediction") 

我該如何爲此添加置信度?

+0

的可能的複製[如何計算斜率的95%置信區間在R中的線性迴歸模型中(http://stackoverflow.com/questions/15180008/how-to-calculate-the-95-confidence-interval-for-the-slope-in​​-a-linear- regressio) – majom

+0

它不是上面提到的重複。 OP是要求一個置信區間而不是置信區間。 – Alex

+0

@Zheyuan Li我也不會認爲它是重複的。首先,術語被使用得很奇怪。如果有人正在尋找它並使用正確的術語,他可能不會找到它。其次,我們在這裏有一個線性迴歸的特例,在其他答案中沒有提到。令人驚訝的是(對我來說),如果你尋找'[r]置信區間線性迴歸',你將不會得到好的結果。所以我認爲這個問題很好。但我也確信在SO上有很多解決方案。但是如果你找不到它們,他們又能怎麼幫忙? – Alex

回答

3

以下代碼將爲您生成好看的迴歸圖。我在代碼中的評論應該清楚地解釋一切。該代碼將使用value,model在您的問題。

## all date you are interested in, 4 years with observations, 10 years for prediction 
all_date <- seq(as.Date("2012-12-31"), by="1 year", length.out = 14) 

## compute confidence bands (for all data) 
pred.c <- predict(model, data.frame(date=all_date), interval="confidence") 

## compute prediction bands (for new data only) 
pred.p <- predict(model, data.frame(date=all_date[5:14]), interval="prediction") 

## set up regression plot (plot nothing here; only set up range, axis) 
ylim <- range(range(pred.c[,-1]), range(pred.p[,-1])) 
plot(1:nrow(pred.c), numeric(nrow(pred.c)), col = "white", ylim = ylim, 
    xaxt = "n", xlab = "Date", ylab = "prediction", 
    main = "Regression Plot") 
axis(1, at = 1:nrow(pred.c), labels = all_date) 

## shade 95%-level confidence region 
polygon(c(1:nrow(pred.c),nrow(pred.c):1), c(pred.c[, 2], rev(pred.c[, 3])), 
     col = "grey", border = NA) 

## plot fitted values/lines 
lines(1:nrow(pred.c), pred.c[, 1], lwd = 2, col = 4) 

## add 95%-level confidence bands 
lines(1:nrow(pred.c), pred.c[, 2], col = 2, lty = 2, lwd = 2) 
lines(1:nrow(pred.c), pred.c[, 3], col = 2, lty = 2, lwd = 2) 

## add 95%-level prediction bands 
lines(4 + 1:nrow(pred.p), pred.p[, 2], col = 3, lty = 3, lwd = 2) 
lines(4 + 1:nrow(pred.p), pred.p[, 3], col = 3, lty = 3, lwd = 2) 

## add original observations on the plot 
points(1:4, rev(value), pch = 20) 

## finally, we add legend 
legend(x = "topleft", legend = c("Obs", "Fitted", "95%-CI", "95%-PI"), 
     pch = c(20, NA, NA, NA), lty = c(NA, 1, 2, 3), col = c(1, 4, 2, 3), 
     text.col = c(1, 4, 2, 3), bty = "n") 

regression plot

的JPEG由代碼生成:

jpeg("regression.jpeg", height = 500, width = 600, quality = 100) 
## the above code 
dev.off() 
## check your working directory for this JPEG 
## use code getwd() to see this director if you don't know 

正如你可以從圖中看到,

  • 置信帶變寬了,當你試圖做預測遠離你觀察的數據;
  • 預測間隔寬於置信區間。

如果您想了解更多關於predict.lm()如何計算內部置信度/預測時間間隔,請閱讀How does predict.lm() compute confidence interval and prediction interval?,並在那裏回答我的答案。

感謝Alex的演示,簡單的使用了visreg包;但我仍然喜歡使用R基地。

+0

我會說這些都沒有信心樂隊。 – Alex

+0

我覺得你用與我不同的方式使用這些術語。對於預測,我將其稱爲預測樂隊,因爲維基百科的確稱爲https://en.wikipedia.org/wiki/Confidence_and_prediction_bands#Prediction_bands。術語置信區間I將用於迴歸參數的區間,如下所述:https://en.wikipedia.org/wiki/Confidence_interval。 – Alex

+0

@ZheyuanLi可以將這兩個地塊作爲一個整體展示,作爲一個連續體?我的意思是,從你的底部情節開始,然後繼續在你的頂部情節? – Oposum

2

您可以簡單地使用visreg::visreg

library(visreg) 
visreg(model) 

enter image description here

如果你有興趣的值:

> head(visreg(model)$fit) 
     date value visregFit visregLwr visregUpr 
1 2012-12-31 13434.5 10753.10 9909.073 11597.13 
2 2013-01-10 13434.5 10807.81 9974.593 11641.02 
3 2013-01-21 13434.5 10862.52 10040.033 11685.00 
4 2013-02-01 13434.5 10917.22 10105.389 11729.06 
5 2013-02-12 13434.5 10971.93 10170.658 11773.21 
6 2013-02-23 13434.5 11026.64 10235.837 11817.44 
+0

謝謝,爲什麼Y軸有5位數值? – Oposum

+0

對不起,我的意思是,爲什麼Y軸與前面的圖相比有不同的值? – Oposum

+0

@ Oposum我認爲這是因爲其他情節使用了未來的數據。但如果你這樣做,你不應該使用一個自信的樂隊,而是一個預測樂隊。置信區間與用於估計迴歸線的數據有關。而預測頻帶應該用於新數據。 – Alex