本文主要是通过交互项，和信息准则AIC,BC以及显著性来确定最优的回归模型：
-------------------------------------
library(MASS)
library(leaps)
## 1: 导入数据：主要变量为E1-E4,剩余变量为：G1-G20
datas<-read.csv(file="data.csv",header=TRUE)
datas
## 2: 展示部分数据集
library(knitr)
kable( head(datas[,1:12]), caption='Partial Data Set' )
## 3: 分析只有主要变量E1-E4来建立关于Y的回归模型
M_E<-lm(Y~E1+E2+E3+E4,data = datas)
summary(M_E)
##### 同时考虑到数据量以及当采用3阶交互时，存在*度过高，无法建立回归模型故采用2阶交互
## 4: 建立在2阶交互项的回归模型
M_raw<-lm(Y~(E1+E2+E3+E4+G1+G2+G3+G4+G5+G6+G7+G8+G9+G10+G11+G12+G13+G14+G15+G16+G17+G18+G19+G20)^2, data=datas )
summary(M_raw)
## 5 此模型包括模型中所有交互项
## 检验残差图
plot(resid(M_raw) ~ fitted(M_raw), main='Residual Plot')
#library(FinTS)
#ArchTest((resid(M_raw))^2,lag=12)
## 6：Box-Cox检验
boxcox(M_raw)
## 7 Box-Cox变换后,建立回归模型
lambda<-0.5
M_trans <- lm(I((Y^lambda-1)/lambda) ~ (.)^2, data=datas )
summary(M_trans)
## 8：新的残差图
plot(resid(M_trans) ~ fitted(M_trans), main='New Residual Plot')
## 9逐步回归
M <- regsubsets( model.matrix(M_trans)[,-1], I((datas$Y^lambda-1)/lambda),
nbest = 1 , nvmax=5,
method = 'forward', intercept = TRUE )
temp <- summary(M)
## 10：软件包为我们提供具有选择性回归模型
Var <- colnames(model.matrix(M_trans))
M_select <- apply(temp$which, 1,
function(x) paste0(Var[x], collapse='+'))
data.frame(cbind( model = M_select, adjR2 = temp$adjr2, BIC = temp$bic))
## 11：确保模型中的重要主效应
M_main <- lm(I((datas$Y^lambda-1)/lambda) ~ ., data=datas)
temp <- summary(M_main)
temp$coefficients[ abs(temp$coefficients[,4]) <= 0.05, ]
# 12 通过使用模型请求中的第二个幂来考虑二阶交互
M_2nd <- lm( I((datas$Y^lambda-1)/lambda) ~ (.)^2, data=datas)
temp <- summary(M_2nd)
temp$coefficients[ abs(temp$coefficients[,4]) <= 0.12, ]
# 13:进一步检验
M_2stage <- lm( I((datas$Y^lambda-1)/lambda) ~ (E2+E3+E4+G12+G13)^2, data=datas)
temp <- summary(M_2stage)
temp$coefficients[ abs(temp$coefficients[,4]) <= 0.05, ]
# 14:最后模型
G12<-datas$G12
G13<-datas$G13
E2<-datas$E2
E3<-datas$E3
E4<-datas$E4
E2_E3<-E2*E3
G12_G13<-G12*G13
Y<-datas$Y
M_2stage <- lm( I((Y^lambda-1)/lambda) ~ E2+E4+E2_E3+G12_G13, data=datas)
summary(M_2stage)