机器学习之线性回归模型预测海藻数量

           线性回归模型使用数据拟合出各个属性之间的函数关系。通过线性回归不仅可以填补缺失值，还可以用来对相关的属性做预测。

#数据集地址：http://www.dcc.fc.up.pt/~ltorgo/DataMiningWithR/datasets2.html
#数据可视化--------------------------------------------------------------------------------------------------
#加载数据两种方式
#第一种：
install.packages(DMwR)
library(DMwR)
head(algae)
#第二种方式
A<-read.table("Analysis.txt",header=F,dec="."
,col.names=c("season","size","speed","mxPH","mn02","C1","NO3","NH4","oPO4","PO4","Chla","a1","a2","a3","a4","a5","a6","a7"),
na.strings=c('xxxxxx'))
#使用mxPH值画直方图
hist(as.numeric(algae$mxPH),prob=T)
#使用qq图看数据是否服从正太分布
install.packages("car")
library(car)

par(mfrow=c(1,2))#将画布分成1行两列的区域
hist(as.numeric(algae$mxPH),prob=T,xlab="",main="Histogram of maximum Ph value",ylim=0:1)#ylim限制y的取值分为0到1
lines(density(as.numeric(algae$mxPH,na.rm=T)))#描绘密度图
rug(jitter(as.numeric(algae$mxPH)))#rug绘图，jitter进行微调避免相同值在视图上的覆盖
qq.plot(as.numeric(algae$mxPH),main="Normal QQ plot of maximum ph")#绘制QQ图
par(mfrow=c(1,1))
#盒图从上往下的线分别含义是（三分位数+1.5*四分位距）、（三分位数）、（均值）、（中位数）、（一分位数）、（一分位数-1.5*四分位距）
boxplot(as.numeric(algae$oPO4),ylab="Orthophosphate(oPO4)")#绘制盒图
rug(jitter(as.numeric(algae$oPO4)),side=2)
abline(h=mean(as.numeric(algae$oPO4,na.rm=T),lty=2))#在均值的位置绘制一条水平线。

#使用盒图观测离群值
plot(as.numeric(algae$NH4),xlab="")
abline(h=mean(as.numeric(algae$NH4,na.rm=T),lty=1))
abline(h=mean(as.numeric(algae$NH4,na.rm=T)+sd(as.numeric(algae$NH4,na.rm=T))),lty=2)
abline(h=median(as.numeric(algae$NH4,na.rm=T)),lty=3)
identity(as.numeric(algae$NH4))
#
plot(as.numeric(algae$NH4),xlab='')
clicked.lines<-identity(as.numeric(algae$NH4))
algae[clicked.lines,]
#
library(lattice)
bwplot(size~a1,data=algae,ylab='River Size',xlab='Algal A1')
#
library(Hmisc)
bwplot(size~a1,data=algae,panel=panel.bpplot,probs=seq(0.01,0.49,by=0.1),datadensity=T,ylab="River Size",xlab="A1")
#
min02<-equal.count(na.omit(algae$mn02),number=4,overlap=1/5)
stripplot(season~a3|min02,data=algae[!is.na(algae$mn02),])
#数据缺失处理方式-------------------------------------------------------------------------------------------------------------
#1.若缺失值很少，直接去除缺失值
#2.取均值、中位数、众数等体现中心趋势的值填补缺失值
#3.通过 变量的相关关系来填补缺失值,最常见的方法就是使用线性回归。步骤：1)先找出相关性较大的变量。2）在相关性较高的变量之间运用显线性回归，建立线性方程求解缺失值
cor(algae[,4:18],use="complete.obs")#产生变量之间的相关值矩阵
symnum(cor(algae[,4:18],use="complete.obs"))#该函数用于改善输出形式
data(algae)
algae <-algae[-manyNAs(algae),]
lm(PO4~oPO4,data=algae)#使用线性回归，建立OP4与oPO4之间的线性方程。OP4=42.897+1.293*oPO4
algae[28,"PO4"]<-42.891+1.293*algae[28,"oOP4"]#使用线性方程填补缺失值
#可以将线性方程写入一个函数中，供以后调用
fillPO4 <- function(oP){
    if(is.na(oP)){
        return(NA)
    }else{
        return(42.897 + 1.293 * oP)
    }
}
algae[is.na(algae$PO4),"PO4"] <- sapply(algae[is.na(algae$PO4),"PO4"],fillPO4)
#绘制条件直方图
histogram(~mxPH|season,data=algae)
algae$season <- factor(algae$season,levels=c("spring","summer","autumn","winter"))
#绘图
histogram(~mxPH|size,data=algae)
histogram(~mxPH|size*speed,data=algae)
stripplot(size~mxPH|speed,data=algae,jitter=T)
#通过探索案例之间的相似性来填补缺失
#使用欧氏距离找出topK个最接近的值，然后用最接近值的均值或者总数替代NA值
#使用高斯核函数设置属性的权重#获取预测模型

data(algae)
algae <-algae[-manyNAs(algae),]
#使用的knn方法中，会建立线性方程填补缺失值。该方法会找出前K个与目标对象相似的对象，
#并用着K个对象的值来填补目标对象的缺失值
clean.algae<-knnImputation(algae,k=10)
#接下来建立用于预测海藻频率的线性回归模型
#第一个参数给出了模型的函数形式，函数的形式使用数据中的其他所有变量来预测变量a1;
#第一个参数中的点“.”代表数据框中的所有除a1外的变量；如果用预测变量mxPH和NH4来预测变量a1,就可以定义模型为a1~mxPH+NH4
#data是R公式al~.模型用到的数据集
#lm的结果是含有线性预测信息的对象
lm.al <-lm(a1~.,data=clean.algae[,1:12])
#使用summary查看lm.al模型信息
#summary中有残差，残差应该是均值为0并且服从正态分布图，残差越小越好
#对于每个多元线性回归方程的系数,R显示他的估计值和标准误差。当然模型中有些参数是不重要的，有些事重要的，因此需要检验系数存在的必要性
#可以进行系数为0的假设检验，即H0：Bi=0.通常使用t检验来验证这些假设，R计算t值，该值为估计系数值与其标准误差的比。
#R将显示与系数相关联的一列Pr(>|t|),表示系数为0这一假设被拒绝的概率
#另一个R输出的诊断信息为R2，R2越接近1表示拟合的越好，否则说明拟合的越差，需要考虑回归模型中的参数的数量
#p值用于表示假设解释变量与目标变量没有依赖关系这一假设的概率
summary(lm.al)
#向后消元法
#使用anova(）来精简线性模型。anova应用的简单线性模型时，该函数提供一个模型拟合的方差序贯分析，也就是说随着公式中
#项数的增加，模型的残差平方和减少
anova(lm.al)
#Pr(>F)大小表示变量对减少模型拟合误差的贡献，该值越大贡献越小
lm2.al <- update(lm.al, . ~ . - season)
#update用于对已有线性模型进行微小的调整，上面用update从lm.al中移除season得到一个新的模型
#模型还是不太理想，使用anova对两个模型进行正式的比较,比较两个模型的相似性Pr(>F)
anova(lm.a1,lm2.a1)
#不断重复上诉过程，剔除属性向后消元，有一种比较简单的向后消元，使用step方法
final.lm<-step(lm.al)
summary(final.lm)
#可以看到R2还是不太理想，因此这种模型是不猛用于al海藻预测分析的
#------------------------------------------------------------------------------------
#R中的另一种回归模型——回归树
library(rpart)
data(algae)
algae<-algae[-manyNAs(algae),]
#建立回归树,当给定的条件满足时构建过程就停止
#1）偏差的减少小于某一个给定界限值时2）当节点中样本数量小于某个给定界限时3）当树的深度大于一个给定的界限值
#上面三个参数分别有rpart函数的三个参数cp、minsplit、maxdepth三个参数指定。他们默认值分别是0.01,20,30
#rpart中实现了一种称为复杂度损失修建的方法，这种修剪法试图估计cp值以确保达到预测的准确性和树的大小之间的最佳折中
rt.al<-rpart(a1 ~ .,data=algae[,1:12])
#显示树
plot(rt.al)
text(rt.al)
prettyTree(rt.al)
#给出一个由函数rpart建立的回归树，R可以生成这些树的一些子树，并估计这些树的性能，可以用printcp方法得到
printcp(rt.al)
#printcp可以看到rpart建议的树，默认树为cp为0.01的树，如果要自己选择一棵树，可以通过使用不同的cp值来建立
rt2.al <- prune(rt.al,cp=0.08)
#rpartXse()函数可以自动运行该过程，他的参数se默认为1
rt.al <- rpartXse(a1 ~ . , data = algae[,1:12])
#使用函数snip.rpart来交互的对树进行检修，有两种方式，第一种是指定需要修建的那个地方的节点号（通过输出树对象来得到节点号
#
first.tree<-rpart(a1 ~. ,data=algae[,1:12])
snip.rpart(first.tree,c(4,7))

#---------------------------------------------------------------------------
#模型的评价与选择——评价模型性能
#回归模型的预测性能是通过将目标变量的预测值与实际值进行比较得到的，并从这些比较重计算平均误差
#一种度量方法是平均绝对误差MAE
#第一步：获取需要评价模型预测性能的预测集个案的预测值，R中要获得模型的预测值，需要使用predict（）函数进行预测
#predict接受两个参数，第一个采纳数是需要应用的模型，另一个参数是数据的测试集，输出结果为相应模型的预测值
lm.predictions.al<-predict(final.lm,clean.algae)
rt.predictions.al<-predict(rt.al,algae)
#第二步：得到预测值之后可以得到平均误差
mae.al.lm<-mean(abs(lm.predictions.al-algae[,"a1"]))
mae.al.rt<-mean(abs(rt.predictions.al-algae[,"a1"]))
#另一种度量误差的度量是均方误差MSE
mes.al.lm<-mean((lm.predictions.al-algae[,"a1"])^2)
mes.al.rt<-mean((rt.predictions.al-algae[,"a1"])^2)
#第三步：使用标准化后的平均绝对误差NMSE，这一统计量是计算模型预测性能和基准模型的预测性能之间的比率，
#通常采用目标变量的平均值来作为基准模型
nmse.al.lm<-mean((lm.predictions.al-algae[,"a1"])^2)/mean((mean(algae[,"a1"])-algae[,"a1"])^2)
nmse.al.rt<-mean((rt.predictions.al-algae[,"a1"])^2)/mean((mean(algae[,"a1"])-algae[,"a1"])^2)

#regr.eval()用来计算线性回归模型的性能指标
regr.eval(algae[,"a1"],rt.predictions.al,train.y=algae[,"a1"])
#两种模型的预测值的可视化分析
old.par<-par(mfrow=c(1,2))
plot(lm.predictions.al,algae[,"a1"],main="Linear Model",xlab="Predictions",ylab="True Values")
abline(0,1,lty=2)
plot(rt.predictions.al,algae[,"a1"],main="Regression Tree","xlab"="Predictions","ylab"="True values")
abline(0,1,lty=2)
par(old.par)
#过滤预测值小于0的数据
sensible.lm.predictions.al<-ifelse(lm.predictions.al < 0,0,lm.predictions.al)
regr.eval(algae[,"a1"],lm.predictions.al,stats=c("mae","mse"))
regr.eval(algae[,"a1"],sensible.lm.predictions.al,stats=c("mae","mse"))

#k折交叉验证——获得模型性能可靠估计的一种常用方法，适用于小数据集
#步骤1：获得K个同样大小的随机训练数据子集，用这k个子集的每一个子集，用出去它之外的其余k-1个子集建立模型，
#然后用第k个子集来评估这个模型，最后存储这个模型的性能指标
#步骤2：对于其余的每个子集重复以上过程，最后又k个性能指标的测量值，这些性能指标是通过在
#没有用于剑魔的数据上计算得到的，这是关键之处
#k折交叉验证估计是这k个性能指标的平均。常见的选择是K=10
#总之遇到一项预测任务时，需要以下决策：
#1.为预测任务选择模型（统一算法的不同参数设定也认为是不同的模型）
#2.选择比较模型性能的评估指标
#3.选择获取评估指标的可靠估计的试验方法


#下面定义用于模型学习和测试的函数,resp函数用于根据公式获得数据集的不变变量值
cv.rpart<-function(form,train,test,...){
  m<-rpartXse(form,train,...)
  p<-predict(m,test)
  mse<-mean((p-resp(form,test))^2)
  c(nmse=mse/mean((mean(resp(form,train)-resp(form,test))^2)))
}

cv.lm<-function(form,train,test,...){
    m<-rpartXse(form,train,...)
    p<-predict(m,test)
    p<-ifelse(p<0,0,p)
    mse<-mean((p-resp(form,test))^2)
    c(nmse=mse/mean((mean(resp(form,train)-resp(form,test))^2)))
}
#定义好用于模型学习和测试的函数后，下面进行交叉验证比较
res<-experimentalComparison(
   c(dataset(a1 ~ .,clean.algae[,1:12],"a1")),
   c(variants("cv.lm"),variants("cv.rpart"),se=c(0,0.5,1)),
   cvSettings(3,10,1234)
)