R語言 線性回歸分析實例 《回歸分析與線性統計模型》page72

y,X1,X2,X3 分別表示第 t 年各項稅收收入(億元)，某國生產總值GDP(億元)，財政支出(億元)和商品零售價格指數(%).

(1) 建立線性模型：

 ① 自己編寫函數：

> library(openxlsx)
> data = read.xlsx("22_data.xlsx",sheet = 1)
> x = data[,-c(1,2)]
> x = cbind(rep(1,17),x)
> x_mat  = as.matrix(x)
> y =matrix(data[,2],ncol = 1)
> res = solve(t(x_mat)%*%x_mat)%*%t(x_mat)%*%y
> res
                    [,1]
rep(1, 17) 19412.8597818
X1             0.2679605
X2            -0.2874013
X3          -297.3653736

 所以各參數的估計值分別為 

 ② lm函數

> lm(y~x_mat)

Call:
lm(formula = y ~ x_mat)

Coefficients:
    (Intercept)  x_matrep(1, 17)          x_matX1  
19412.859781545               NA      0.267960511  
        x_matX2          x_matX3  
   -0.287401287   -297.365373557  

 於是各參數的估計值分別為

這兩個方法的結果是一樣的。

 

(2)要求實驗報告中畫出矩陣散點圖，給出參數的點估計、區間估計、t檢驗值、判定系數和模型F檢驗的方差分析表

繪制矩陣散點圖。

library(graphics)
pairs(data[,-1]pch = 21,bg = c('red','green3','blue'))
# pch參數是控制點的形狀，bg是控制點的顏色

 

下面代碼給出參數的點估計，t檢驗值，判定系數

> summary(lm(y~x_mat+1))

Call:
lm(formula = y ~ x_mat + 1)  #調用

Residuals:      #殘差統計量，殘差第一四分位數(1Q)和第三分位數(3Q)有大約相同的幅度，意味着有較對稱的鍾形分布
    Min      1Q  Median      3Q     Max  
-4397.9 -1102.4   153.8  1184.4  2934.6 

Coefficients: (1 not defined because of singularities)        
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)      1.941e+04  3.524e+04   0.551    0.591    
x_matrep(1, 17)         NA         NA      NA       NA    
x_matX1          2.680e-01  4.466e-02   6.000 4.45e-05 ***
x_matX2         -2.874e-01  1.668e-01  -1.723    0.109    
x_matX3         -2.974e+02  3.688e+02  -0.806    0.435    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#標記為Estimate的列包含由最小二乘法計算出來的估計回歸系數。
#標記為Std.Error的列是估計的回歸系數的標准誤差。
#從理論上說，如果一個變量的系數是0，那么該變量將毫無貢獻。然而，這里顯示的系數只是估計，它們不會正好為0.
#因此，我們不禁會問：從統計的角度而言，真正的系數為0的可能性有多大？這是t統計量和P值的目的，在匯總中被標記為t value和Pr(>|t|)
#P值估計系數不顯著的可能性，有較大P值的變量是可以從模型中移除的候選變量

Residual standard error: 2013 on 13 degrees of freedom
Multiple R-squared:  0.9982,	Adjusted R-squared:  0.9977 
F-statistic:  2348 on 3 and 13 DF,  p-value: < 2.2e-16

#Residual standard error 表示殘差的標准差，F-statistic 表示F的統計量

　　區間估計？方差分析表？

#參數的區間估計
fm = lm(y~x_mat)
confint(fm,level = 0.95)
anova(fm) #方差分析

 

(3)保留模型中線性關系顯著的預測變量確定最后的模型，並利用R軟件中的"predict"語句預測2017年的稅收收入

 根據回歸分析結果，只有變量X1具有顯著性。所以模型中僅保留變量X1。

構造模型：

x_mat = cbind(rep(1,17),data[,3]) 
y = data[,2]
res = lm(y~x_mat)
res
> res

Call:
lm(formula = y ~ x_mat)

Coefficients:
(Intercept)       x_mat1       x_mat2  
 -6213.0189           NA       0.1915 

 

 該模型為：Y = -6213.0189 + 0.1915 X1

接下來預測2017年的稅收收入，先根據數據data對 t 和 y 之間的關系進行回歸分析

t = data[,1]
y = data[,2]
res = lm(y~t)
res

> res

Call:
lm(formula = y ~ t)

Coefficients:
(Intercept)            t  
  -16428607         8213  

 所以 t 與 y 的關系為：y  = -16428607 + 8213 t

 預測 2017 年的稅收收入：

> newdata = data.frame(t = 2017)
> pre = predict(res,newdata,interval = "prediction",level = 0.95)
> pre
       fit      lwr      upr
1 136337.8 116018.1 156657.4

　　

（4）方差齊性檢驗，正態性檢驗，誤差相關性的DW檢驗

rm(list = ls())
library(openxlsx)
data = read.xlsx("22_data.xlsx",sheet = 3)
data
attach(data) #執行此命令之后，直接用列名引用數據
fm = lm(y~X1+X2+X3)
summary(fm)

　　

#殘差圖：方差齊性檢驗
ei = resid(fm) #殘差
X = cbind(1,as.matrix(data[,3:5])) #創建設計陣，注意as.matrix的對象
t = ti(ei,X) #外部學生化殘差
r = ri(ei,X) #學生化殘差
plot(fitted(fm),t,xlab = "y估計值",ylab = "學生化殘差",ylim = c(-4,4))
abline(h = 0)

　　

1)上圖表示，圍繞x軸上下波動，不是均勻分布，方差不齊，但是這個殘差圖代表什么呢？
2)第17個點是異常點

#正態概率圖：正態性檢驗
plot_ZP(t)

　　

 

 

 除了17號點，其他的點基本在一條直線。

plot_ZP(t[1:16])

　　

 

 

#誤差相關性的DW檢驗
library(lmtest)
dw = dwtest(fm) #DW檢驗函數

> dw

Durbin-Watson test

data: fm
DW = 0.90086, p-value = 0.0006763
alternative hypothesis: true autocorrelation is greater than 0

p值遠小於0.05，拒絕原假設ρ=0，所以誤差是自相關的

plot(ei) #繪制時間序列中的殘差圖
abline(h = 0)

　　

圖像表明：誤差是正自相關的。
正自相關：如果前一個殘差大於0，那么后一個殘差大於0的概率較大
負自相關：富國前一個殘差大於0，那么后一個殘差小於0的概率較大

（4）強影響點分析，異常點分析

#強影響點分析
#方法：手指律、刀切法、cook距離、deffits
#influence.measures(fm)
influence.measures(fm)
# 函數得到的回歸診斷共有9列，
# 其中1,2,3,4列是dfbetas值（對應常數與變量x），
# 第5列是dffits的准則值，
# 第6列是covratio的准則值,
# 第7列是cook值,第8列是帽子值（高杠桿值），
# 第9列是影響點的標記，
# inf表明16，17號是強影響點。

　　

#cook距離判斷強影響點
res = cooks.distance(fm)
> res[res>4/(17-3-1)]
      17 
1.266196 
> #17號是強影響點

　　

#異常點
#方法：dfbetas、F統計量、outlierTest()
library(car)
outlierTest(fm)

    rstudent unadjusted p-value Bonferroni p
17 -4.277398          0.0010739     0.018257

Bonferroni p < 0.05 , 結果顯示17號點是異常點

　　

#使用influencePlot()將異常點繪入圖中
influencePlot(fm)

　　

 

 

#F統計量page88
calcu_F = function(p,r) #p回歸參數個數，r學生化殘差
{
  n = length(r)
  ans = (n-p-2)*(r**2)/(n-p-1-r**2)
  return(ans)
}
ff = calcu_F(3,r)
#與自由度為1,n-p-2,顯著性水平為α的F值比較
df_val = qf(1-0.05,1,12)

ff[ff>df_val] #檢驗顯著性
      17 
18.29613 

　　

（6）模型失擬檢測

模型失擬檢測
1.有重復值用失擬檢測的正規檢驗《線性回歸導論》
2.無重復值
1）當模型的預測變量個數多余1時，考慮偏殘差法
2）無論預測變量個數，近鄰點估計純誤差方法都可以

# 偏殘差：
data
beta = coef(fm)
beta1 = beta[2]
beta2 = beta[3]
beta3 = beta[4]
#第一個預測變量：
par(mfrow = c(2,2))
cancha1 = ei + beta1*(data$X1)
plot(data$X1,cancha1)
#第二個預測變量：
cancha2 = ei + beta2*(data$X2)
plot(data$X2,cancha2)
#第三個預測變量：
cancha3 = ei + beta3*(data$X3)
plot(data$X3,cancha3)
#通過圖像，第三個預測變量不🆗不成線性關系
detach(data)