統計學習的模型與操作一、預測變量為連續型數值變量讀入數據：rm(list=ls())setwd作工作\\課程\\經管軟件應用\2016")install.packages("mboost")library("mboost")data("bodyfat",package="TH.data")str(bodyfat)write.table(bodyfat,file="bodyfat.csv",s=FALSE)數據結構及描述性統計Str(bodyfat)Summary(bodyfat)將數據集分為訓練集和測試集set.seed(1234)ind=sample(2,nrow(bodyfat),replace=TRUE,prob=c(0.7,0.3))bodyfat_train=bodyfat[ind==1,]bodyfat_test=bodyfat[ind==2,]對測試集進行廣義線性模性預測###構建廣義線性模性myformula=DEXfat~age+waistcirc+hipcirc+elbowbreadth+kneebreadthbodyfat.glm=glm(myformula,data=bodyfat_train,family=gaussian("log"))summary(bodyfat.glm)###繪圖，比較觀測值與預測值library(ggplot2)###訓練集中的觀測值與預測值。p=qplot(bodyfat_train$DEXfat,pred1,colour="red",geom=c("point"))p=p+geom_abline(intercept=0,slope=1,colour="blue")p###測試集中的觀測值與預測值pred2=predict(bodyfat.glm,newdata=bodyfat_test,type="response")p=qplot(bodyfat_test$DEXfat,pred2,colour=I("green"),geom=c("point"))p=p+geom_abline(intercept=0,slope=1,colour="blue")+ggtitle(‘測試集的觀測值與預測值(廣義線性模型)")p###測試集中的預測總誤差。error1=sum((pred2-bodyfat_test$DEXfat)A2)###測試集中的均方根誤差RMS1=sqrt(mean((pred2-bodyfat_test$DEXfat)A2))對測試集進行決策樹分析###構建決策樹模型library(rpart.plot)library(rpart)bodyfat.rpart=rpart(myformula,data=bodyfat_train,control=rpart.control(minsplit=10))summary(bodyfat.rpart)###繪制決策樹模型圖plot(bodyfat.rpart,uniform=TRUE)text(bodyfat.rpart,use.n=T)prp(bodyfat.rpart)###測試集中的觀測值與預測值pred3=predict(bodyfat.rpart,newdata=bodyfat_test)p=qplot(bodyfat_test$DEXfat,pred2,colour=I("green"),geom=c("point"))p=p+geom_abline(intercept=0,slope=1,colour="blue")+ggtitle(‘測試集的觀測值與預測值(決策樹模型)")p###模型預測能力評價error2=sum((pred3-bodyfat_test$DEXfat)A2)RMS2=sqrt(mean((pred3-bodyfat_test$DEXfat)A2))對預測集進行隨機森林分析###隨機森林library(randomForest)bodyfat.ren=randomForest(myformula,data=bodyfat_train)summary(bodyfat.ren)pred4=predict(bodyfat.ren,newdata=bodyfat_test)error3=sum((pred4-bodyfat_test$DEXfat)A2)RMS3=sqrt(mean((pred4-bodyfat_test$DEXfat)A2))p=qplot(bodyfat_test$DEXfat,pred4,colour=I("red"),geom=c("point"))p=p+geom_abline(intercept=0,slope=1,colour="blue")+ggtitle(‘測試集的觀測值與預測值(隨機森林模型)")p二、預測變量為二分式分類變量rm(list=ls())setwd("E:\\E\\yang\\工作\\課程\\經管軟件應用\\2016")#Readinthedatacensus=read.csv("census.csv")str(census)head(census)#splitdata(兩個函數，caTools包中的sample.split函數，或者base包中的sample函數)library(caTools)set.seed(2000)spl=sample.split(census$over50k,SplitRatio=0.6)train=subset(census,spl==TRUE)test=subset(census,spl==FALSE)set.seed(2000)ind=sample(2,nrow(census),replace=TRUE,prob=c(0.6,0.4))train1=census[ind==1,]summary(train)summary(train1)#logit(logit模型預測和分析)trainlogit=glm(over50k~.,data=train,family="binomial")summary(trainlogit)#predictionlogitpre=predict(trainlogit,newdata=test,type="response")A1=table(test$over50k,logitpre>0.5)sum(diag(A1))/nrow(test)table(test$over50k)#auclibrary(ROCR)ROCRpredTest=prediction(logitpre,test$over50k)perf1=performance(ROCRpredTest,"tpr","fpr")plot(perf1)auc=as.numeric(performance(ROCRpredTest,"auc")@y.values)auc#treelibrary(rpart)library(rpart.plot)CARTcensus=rpart(over50k~.,data=train,method="class")prp(CARTcensus)cartpre=predict(CARTcensus,newdata=test,type="class")table(test$over50k,cartpre)predictTest=predict(CARTcensus,newdata=test)predictTest2=predictTest[,2]#ComputetheAUC:ROCRpred=prediction(predictTest2,test$over50k)perf2=performance(ROCRpred,"tpr","fpr")plot(perf2)as.numeric(performance(ROCRpred,"auc")@y.values)head(predictTest)A2=table(test$over50k,predictTest2>0.5)sum(diag(A2))/nrow(test)#RANDOMFORESTlibrary(randomForest)Forest=randomForest(over50k~.,data=train)forestpre=predict(Forest,newdata=test)A3=table(test$over50k,forestpre)sum(diag(A3))/nrow(test)#R