利用 R 进行多元回归估计
楚新元 / 2021-08-23
案例来源于古扎拉蒂《经济计量学精要》(第四版)第二章例 2-5 古董钟与拍卖价格。之所以选择此案例是因为此案例只有 3 个变量，32 个样本，便于手工验证。
读取数据
library(readxl)
library(dplyr)
library(kableExtra)
data = read_excel("data/Table2_14.xls", skip = 5)
data %>%
kable(
caption = "古董钟与拍卖价格",
format = "html"
) %>%
kable_styling(
bootstrap_options = "striped",
fixed_thead = TRUE,
font_size = 14,
) %>%
column_spec(1, width = "80px") %>%
column_spec(2:4, width = "120px")
| Observation | Price | Age | Number of Bidders |
|---|---|---|---|
| 1 | 1235 | 127 | 13 |
| 2 | 1080 | 115 | 12 |
| 3 | 845 | 127 | 7 |
| 4 | 1552 | 150 | 9 |
| 5 | 1047 | 156 | 6 |
| 6 | 1979 | 182 | 11 |
| 7 | 1822 | 156 | 12 |
| 8 | 1253 | 132 | 10 |
| 9 | 1297 | 137 | 9 |
| 10 | 946 | 113 | 9 |
| 11 | 1713 | 137 | 15 |
| 12 | 1024 | 117 | 11 |
| 13 | 2131 | 170 | 14 |
| 14 | 1550 | 182 | 8 |
| 15 | 1884 | 162 | 11 |
| 16 | 2041 | 184 | 10 |
| 17 | 854 | 143 | 6 |
| 18 | 1483 | 159 | 9 |
| 19 | 1055 | 108 | 14 |
| 20 | 1545 | 175 | 8 |
| 21 | 729 | 108 | 6 |
| 22 | 1792 | 179 | 9 |
| 23 | 1175 | 111 | 15 |
| 24 | 1593 | 187 | 8 |
| 25 | 1147 | 137 | 8 |
| 26 | 1092 | 153 | 6 |
| 27 | 1152 | 117 | 13 |
| 28 | 1336 | 126 | 10 |
| 29 | 785 | 111 | 7 |
| 30 | 744 | 115 | 7 |
| 31 | 1356 | 194 | 5 |
| 32 | 1262 | 168 | 7 |
变量含义
- Price：古董钟中标的价格
- Age：钟表的年代
- Number of Bidders：投标人的个数
利用 R 自带函数进行参数估计
fit = lm(Price ~ Age + `Number of Bidders`, data = data)
summary(fit)
##
## Call:
## lm(formula = Price ~ Age + `Number of Bidders`, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -208.60 -119.04 15.73 101.84 212.54
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1336.0493 175.2725 -7.623 2.10e-08 ***
## Age 12.7414 0.9124 13.965 2.09e-14 ***
## `Number of Bidders` 85.7641 8.8020 9.744 1.19e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 134.6 on 29 degrees of freedom
## Multiple R-squared: 0.8906, Adjusted R-squared: 0.8831
## F-statistic: 118.1 on 2 and 29 DF, p-value: 1.161e-14
手动计算参数值和各个统计量(含截距项)
生成矩阵 \(X\)
data %>%
mutate(X1 = 1) %>%
select(X1, Age, `Number of Bidders`) %>%
as.matrix() -> X
print(X)
## X1 Age Number of Bidders
## [1,] 1 127 13
## [2,] 1 115 12
## [3,] 1 127 7
## [4,] 1 150 9
## [5,] 1 156 6
## [6,] 1 182 11
## [7,] 1 156 12
## [8,] 1 132 10
## [9,] 1 137 9
## [10,] 1 113 9
## [11,] 1 137 15
## [12,] 1 117 11
## [13,] 1 170 14
## [14,] 1 182 8
## [15,] 1 162 11
## [16,] 1 184 10
## [17,] 1 143 6
## [18,] 1 159 9
## [19,] 1 108 14
## [20,] 1 175 8
## [21,] 1 108 6
## [22,] 1 179 9
## [23,] 1 111 15
## [24,] 1 187 8
## [25,] 1 137 8
## [26,] 1 153 6
## [27,] 1 117 13
## [28,] 1 126 10
## [29,] 1 111 7
## [30,] 1 115 7
## [31,] 1 194 5
## [32,] 1 168 7
生成矩阵 \(Y\)
Y = as.matrix(data$Price)
print(Y)
## [,1]
## [1,] 1235
## [2,] 1080
## [3,] 845
## [4,] 1552
## [5,] 1047
## [6,] 1979
## [7,] 1822
## [8,] 1253
## [9,] 1297
## [10,] 946
## [11,] 1713
## [12,] 1024
## [13,] 2131
## [14,] 1550
## [15,] 1884
## [16,] 2041
## [17,] 854
## [18,] 1483
## [19,] 1055
## [20,] 1545
## [21,] 729
## [22,] 1792
## [23,] 1175
## [24,] 1593
## [25,] 1147
## [26,] 1092
## [27,] 1152
## [28,] 1336
## [29,] 785
## [30,] 744
## [31,] 1356
## [32,] 1262
计算样本容量 \(n\) 和自由度 \(df\)
n = nrow(data)
k = ncol(X)
df = n - k
print(c(n = n, k = k, df = df))
## n k df
## 32 3 29
估计参数 \(\beta\)
Beta_hat = solve(t(X) %*% X) %*% t(X) %*% Y
colnames(Beta_hat) = "Beta_hat"
print(Beta_hat)
## Beta_hat
## X1 -1336.04929
## Age 12.74138
## Number of Bidders 85.76407
计算回归残差 \(e\)
residuals = Y - (X %*% Beta_hat)
colnames(residuals) = "residuals"
print(residuals)
## residuals
## [1,] -162.03929
## [2,] -78.37862
## [3,] -37.45489
## [4,] 204.96516
## [5,] -119.19094
## [6,] 52.71275
## [7,] 141.22465
## [8,] 49.54599
## [9,] 115.60314
## [10,] 70.39635
## [11,] 17.01874
## [12,] -74.09732
## [13,] 100.31715
## [14,] -118.99505
## [15,] 212.54042
## [16,] 174.99405
## [17,] -146.55296
## [18,] 21.29270
## [19,] -185.71707
## [20,] -34.80536
## [21,] 174.39547
## [22,] 75.46503
## [23,] -189.70529
## [24,] -139.70197
## [25,] 51.36721
## [26,] -35.96679
## [27,] -117.62546
## [28,] 208.99429
## [29,] 106.40725
## [30,] 14.44171
## [31,] -208.59945
## [32,] -142.85161
计算回归标准误 \(\hat{\sigma}\)
Sigma_hat = sqrt(sum(residuals^2) / df)
c(Sigma_hat = Sigma_hat)
## Sigma_hat
## 134.6083
计算 \(\beta\) 估计量的方差与标准误
Beta_hat_var = solve(t(X) %*% X) * (Sigma_hat)^2
Beta1_hat_se = sqrt(Beta_hat_var[1, 1])
Beta2_hat_se = sqrt(Beta_hat_var[2, 2])
Beta3_hat_se = sqrt(Beta_hat_var[3, 3])
c(
Beta1_hat_se = Beta1_hat_se,
Beta2_hat_se = Beta2_hat_se,
Beta3_hat_se = Beta3_hat_se
)
## Beta1_hat_se Beta2_hat_se Beta3_hat_se
## 175.2724965 0.9123559 8.8019949
计算 \(t\) 统计量
t_Beta1_hat = Beta_hat[1, 1] / Beta1_hat_se
t_Beta2_hat = Beta_hat[2, 1] / Beta2_hat_se
t_Beta3_hat = Beta_hat[3, 1] / Beta3_hat_se
c(
t_Beta1_hat = t_Beta1_hat,
t_Beta2_hat = t_Beta2_hat,
t_Beta3_hat = t_Beta3_hat
)
## t_Beta1_hat t_Beta2_hat t_Beta3_hat
## -7.622698 13.965366 9.743708
计算拟合优度 \(R^2\)
TSS = sum((Y - mean(Y))^2)
RSS = sum(residuals^2)
ESS = TSS - RSS
R_squared = 1 - RSS/TSS
c(R_squared = R_squared)
## R_squared
## 0.8906143
计算修正的拟合优度 \(Adjusted\ R^2\)
Adjusted_R_squared = 1 - (1 - R_squared) * (n - 1) / (n - k)
c(Adjusted_R_squared = Adjusted_R_squared)
## Adjusted_R_squared
## 0.8830705
计算 \(F\) 统计量
F_statistic = (R_squared / (k - 1)) /
((1 - R_squared) / (n - k))
c(F_statistic = F_statistic)
## F_statistic
## 118.0585