Hồi qui đa biến
Y=c(12.1, 11.9, 10.2, 8.0, 7.7, 5.3, 7.9, 7.8, 5.5, 2.6)
X1=c(0:9)
X2= c(7, 4, 4, 6, 4, 2, 1, 1, 1, 0)
df=data.frame(Y,X1,X2)
yx1m=lm(Y ~ X1, data=df)
yx2m=lm(Y ~ X2, data=df)
yx1x2m=lm(Y ~ X1 +X2, data=df)
summary(yx1m)
##
## Call:
## lm(formula = Y ~ X1, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1606 -1.0735 0.1742 0.8621 2.0970
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.8545 0.8283 14.312 5.54e-07 ***
## X1 -0.8788 0.1552 -5.664 0.000474 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.409 on 8 degrees of freedom
## Multiple R-squared: 0.8004, Adjusted R-squared: 0.7755
## F-statistic: 32.08 on 1 and 8 DF, p-value: 0.0004737
summary(yx2m)
##
## Call:
## lm(formula = Y ~ X2, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.702 -1.533 -0.034 1.667 3.066
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.0980 1.1222 4.543 0.00189 **
## X2 0.9340 0.2999 3.114 0.01436 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.121 on 8 degrees of freedom
## Multiple R-squared: 0.548, Adjusted R-squared: 0.4915
## F-statistic: 9.698 on 1 and 8 DF, p-value: 0.01436
summary(yx1x2m)
##
## Call:
## lm(formula = Y ~ X1 + X2, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.46078 -0.33384 0.00026 0.81856 1.98476
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.7076 2.9785 4.938 0.00168 **
## X1 -1.2042 0.3614 -3.332 0.01255 *
## X2 -0.4629 0.4642 -0.997 0.35187
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.41 on 7 degrees of freedom
## Multiple R-squared: 0.8252, Adjusted R-squared: 0.7753
## F-statistic: 16.53 on 2 and 7 DF, p-value: 0.002232
AIC(yx1m)
## [1] 39.00848
AIC(yx2m)
## [1] 47.18316
AIC(yx1x2m)
## [1] 39.6801
Nhận xét: mô hình yx1m có AIC=39.00848 nhỏ nhất: mô hình tối
ưu.
Phân tích khả năng nhận dạng mô hình theo phương pháp BMA
library(BMA)
## Loading required package: survival
## Loading required package: leaps
## Loading required package: robustbase
##
## Attaching package: 'robustbase'
## The following object is masked from 'package:survival':
##
## heart
## Loading required package: inline
## Loading required package: rrcov
## Scalable Robust Estimators with High Breakdown Point (version 1.7-7)
# Ma trận biến độc lập và biến phụ thuộc
#X <- df[, c("X1", "X2")]
#y <- df$Y
#bma_model <- bicreg(X, y, strict = FALSE)
bma = bicreg(df[, c("X1", "X2")], df$Y, strict=FALSE,OR=20)
summary(bma)
##
## Call:
## bicreg(x = df[, c("X1", "X2")], y = df$Y, strict = FALSE, OR = 20)
##
##
## 2 models were selected
## Best 2 models (cumulative posterior probability = 1 ):
##
## p!=0 EV SD model 1 model 2
## Intercept 100.0 12.9403 2.3917 11.8545 14.7076
## X1 100.0 -1.0026 0.2993 -0.8788 -1.2042
## X2 38.1 -0.1762 0.3640 . -0.4629
##
## nVar 1 2
## r2 0.800 0.825
## BIC -13.8120 -12.8378
## post prob 0.619 0.381
Nhậ xét:
p!=0: X1: 100% ⇒ luôn xuất hiện trong mô hình tốt nhất.
X2: 38.1% ⇒ chỉ xuất hiện trong một phần mô hình.
R square: Model 2 chỉ tăng nhẹ R² (0.825 so với 0.800) → cải thiện
rất ít,nhưng BIC lại kém hơn ### (ít âm hơn = mô hình kém hơn một
chút).
Vì thế, BMA cho rằng biến X2 không đóng góp nhiều ý nghĩa thống kê,
chỉ có ~38% xác suất nên được ### giữ lại.
Kết luận: Mô hình chính: Y = 11.85 − 0.88·X1
Khi áp dụng Bayesian Model Averaging, ta thấy rằng biến X1 có ảnh
hưởng rõ rệt và ổn định đến Y.
Biến X2 chỉ có ảnh hưởng nhẹ, không nhất quán giữa các mô hình.
