#正式测试（根据自己的数据进行分析）
#下载xlsx函数包
install.packages("xlsx")
#导入xlsx包
library(xlsx)
#测试数据一
#从xlsx表中读取数据
test.data<-read.xlsx("2012_2场1.xlsx",1,header=TRUE,as.data.frame=TRUE)
#打出散点图
plot(test.data)
lm.test<-lm(formula = test.data$density~.,data=test.data)
summary(lm.test)
#我们可以看到除了日期以外，其他各解释变量的P值都比较大，表现为不相关，
#但是可以看到R^2为0.7935还是比较高的，比较理想

#查看各解释变量与被解释变量之间的相关系数
cor(test$density,test$Mtemp)
cor(test$density,test$Htemp)
cor(test$density,test$Ltemp)
cor(test$density,test$huminity)
cor(test$density,test$sunshine)
cor.test(test$density,test$Mtemp)
cor.test(test$density,test$Htemp)
cor.test(test$density,test$Ltemp)
cor.test(test$density,test$huminity)
cor.test(test$density,test$sunshine)
#会发现计算得到的相关系数并不低，只有相对湿度与虫口密度的相关系数为负数
#可以知道相对湿度与虫口密度之间是呈负相关的，是合理的

#多因素方差分析的R函数和示例
result<-aov(test$density~test$Mtemp+test$Htemp+test$Ltemp+test$huminity+test$sunshine)
anova(result)
#通过结论我们可以看出最高温度、平均温度、相对湿度、光照时长都比较理想
#若显著性水平α为0.05，则因为Mtemp、Htemp、huminity、sunshine的概率P小于α
#则应拒绝原假设，认为这几个因素对虫口密度有显著影响
#因为Ltemp的概率P值大于α，所以不能拒绝原假设，认为Ltemp对虫口密度没有显著的影响

#因为实验证明Mtemp、Htemp、huminity、sunshine和虫口密度有一定的相关关系，但是模型并不好
#那么我们现在就来优化模型
#利用残差图判断误差项是否满足高斯—马尔科夫假定
#对于误差项是否满足高斯-马尔科夫假定，可以从三个方面考察：独立性、正态性、等方差性
lm.test<-lm(formula = test.data$density~.,data=test.data)
par(mfrow=c(2,2))
plot(lm.test)

#第一幅图（横坐标为被解释变量的拟合值，纵坐标为残差）等方差性的本质上是指误差项的方差
#不随各解释变量的变化而变化，多元中考察解释变量全体最简洁最有效的方式就是考察他们的线性组合
#也即考察被解释变量的拟合值。如果等方差时，残差的方差不应随拟合值的增大而呈明显的趋势规律性
#在等方差无法满足的情况下，一般的解决方法有三种：
#第一，扔采用原有模型，但利用加权最小二乘估计法估计参数（P217）
#第二，采用异方差——稳健t统计量进行回归系数的显著性检验
#第三，对被解释变量做BOX-COX变换后在建立模型（变换的最常见的形式是取自认对数）
#QQ图是关于残差项的Q-Q图，因仍有部分残差点有规律的落在对角线之外，故可近似认为残差项的
#正态性并未得到满足。这种情况下，通常可先对被解释变量做调整变换后再建立模型

#检验误差项的独立性假定(程辑包‘car’是用R版本3.2.4 来建造的 )
lm.test<-lm(formula = test.data$density~.,data=test.data)
#下载car函数包
install.packages("car")
#导入car函数包
library("car")
durbinWatsonTest(lm.test)

#回归值与残差的残差图
coef(lm.test)

#回归诊断：诊断数据中的异常观测点（探测高杠杆值点、探测离群点、探测强影像点）
#探测高杠杆值点——指的是那些在解释变量方向上取值异常、在被解释变量方向上取值正常的点
#所谓在解释变量上取值异常，是指该观测点在x上的取值远远偏离绝大部分
#观测的平均水平。杠杆值就是对这个远离程度的度量，如果观测的杠杆值大于
#2倍或是3倍的杠杆值得平均值，则可认为该观测值为高杠杆值点
#探测离群点——是指那些在被解释变量方向上取值异常的观测点。离群点在y上的取值远远偏离绝大
#部分观测的取值范围。通常表现为：
#第一，离群点的被解释变量取值非常大或非常小。因被解释变量服从正态分布，而
#这些观测点的被解释变量往往在2个或3个标准差之外取值
#第二，离群点的残差的绝对值非常大。在误差项满足高斯-马尔科夫假定的条件下
#（零均值等方差），离群值的残差一般在2个或3个残差标准差之外取值
#离群点会将回归经验方程“拉向”自己，对回归模型的参数估计造成影响
#探测影像点——强影响点是这样的观测点，包含或剔除该观测点，会导致回归经验方程的截距或
#斜率发生较大的变化。强影响点不一定在被解释变量的方向上取异常值，也不一定
#在解释变量方向上取异常值，但会在两者的比值上 出现异常。
#库克距离是一种探测强影响点的度量方法。
LeveragePlot<-function(fit)
{
　　Np<-length(coefficients(fit))-1 #获得解释变量的个数
　　N<-length(fitted(fit)) #获得样本量
　　plot(hatvalues(fit),main="观测点的杠杆值序列图",ylab="杠杆值",xlab="观测编号")
　　#hatvalues()红帽值和回归删除诊断
　　abline(2*(Np+1)/N,0,col="red",lty=2)
　　abline(3*(Np+1)/N,0,col="red",lty=2)
　　identify(1:N,hatvalues(fit),names(hatvalues(fit)))
}
lm.test<-lm(formula = test.data$density~.,data=test.data)
LeveragePlot(lm.test)
#我们可以有一个高杠杆值点，它比平均杠杆值要高

#探测离群点的统计检验方法是以学生化残差为检验统计量所进行的检验（P225）
#R访问学生化残差的函数是restudent
lm.test<-lm(formula = test.data$density~.,data=test.data)
library("car")
outlierTest(lm.test)
#具体情况参照《基于R的统计分析与数据挖掘》（P226）

#探测强影响点
lm.test<-lm(formula = test.data$density~.,data=test.data)
#获得库克距离
plot(cooks.distance(lm.test),main="Cook's distance",cex=0.5)
Np<-length(coefficients(lm.test))-1
N<-length(fitted(lm.test))
CutLevel<-4/(N-Np-1)
plot(lm.test,which=4)
abline(CutLevel,0,lty=2,col="red")

#异常观测点的综合可视化
install.packages("car")
library("car")
lm.test<-lm(formula = test.data$density~.,data=test.data)
influencePlot(lm.test)

#回归诊断：多重共线性的诊断
#多重共线性的测度统计量：检验多重共线性的统计量有容忍度和方差膨胀因子，互为倒数关系
#容忍度是测度解释变量间多重共线性的重要统计量
#容忍度的取值范围为0~1，越接近1表示多重共线性越弱，越接近0表示多重共线性越强
#所以，线性回归分析中，各个解释变量的容忍度不应太小
lm.test<-lm(formula = test.data$density~.,data=test.data)
library("car")
vif(lm.test)
#出现较严重的多重共线性时的一般处理策略是：剔除其中的某个或某几个解释变量，重新
#建立线性回归模型

#回归建模策略
#多重共线性给出的启示：并非所有与被解释变量有线性关系的解释变量都应该考虑进线性回归模型
#线性回归模型包含多个解释变量，不同的建模策略可能导致最终留在模型中的解释变量不完全相同
#进而最终得到的模型也有差别。选择建模策略可以从以下两个角度考虑：
#第一，所建立的模型应是拟合优度最高的模型，所建立的模型应最大限度的拟合样本数据
#第二，所建立的模型应是带惩罚的拟合优度最高的模型。
#在上述建模策略下还需进一步考虑线性回归的具体实施方案：
#一方面，规定解释变量进出模型的顺序。为避免严重的多重共线性，确定解释变量进出模型的“门槛”
#即评判标准。这取决于解释变量的筛选策略
#另一方面，如何确保在解释变量筛选的框架下，不遗漏可能的最佳模型，即全子集回归问题

#利用AIC（赤池信息量准则）（Akaike information criterion）来进行模型优化
#BIC（贝叶斯信息准则）（Bayesian information criterion）来进行模型优化
lm.test1<-lm(formula = test.data$density~.,data=test.data)
lm.test2<-lm(formula=test.data$density~test.data$Htemper+test.data$Mtemper+test.data$humidity+test.data$sunshine,data=test.data)
#计算赤池信息量准则
AIC(lm.test1,lm.test2)
#计算贝叶斯信息准则
BIC(lm.test1,lm.test2)

#解释变量的筛选策略
#向前引入法：从一元回归开始，逐步增加变量，使指标值达到最优为止
#向后剔除法：从全变量回归方程开始，逐步删去某个变量，使指标值达到最优为止
#逐步筛选法：综合上述两种方法
install.packages("MASS")
library("MASS")
lm.test<-lm(formula = test.data$density~test.data$Ltemp+test.data$Htemper+test.data$Mtemper+test.data$humidity+test.data$sunshine,data=test.data)
summary(lm.test)
stepAIC(lm.test,direction = "backward")#采用赤池信息准则，使用向后剔除法进行优化模型

#优化剔除后的模型
lm.test<-lm(test.data$density ~ test.data$Htemper + test.data$humidity,data=test.data)
summary(lm.test)

#对Y进行开方
lm.test<-lm(sqrt(test.data$density) ~ test.data$Htemper + test.data$humidity,data=test.data)
summary(lm.test)

#对Y进行取对数
lm.test<-lm(log(test.data$density) ~ test.data$Htemper + test.data$humidity,data=test.data)
summary(lm.test)

#利用12、13、14、15年的2场数据总体进行测验
test<-read.xlsx("2012345_2场1.xlsx",1,header=TRUE,as.data.frame=TRUE)#从xlsx表中读取数据
lm.test<-lm(formula = test$density ~ test$Ltemper + test$humidity + test$sunshine, data = test.data)
summary(lm.test)