下面我們開始運用R語言分析來源於UCI數據庫中的關於白酒品質研究的數據集進行算法演示，該數據集是關於白酒中的各項變量對白酒品質的影響情況。

這裡將利用該數據集建立出適合的單隱藏層前饋人工神經網絡模型，並對所建立的模型進行相應的分析，查看建立模型的預測能力如何。

數據探索

我們先從UCI數據庫中下載關於白酒品質的的"winequality-white.csv"數據集，下載地址為：http://archive.ics.uci.edu/ml/machine-learning-databases/wine-quality/。這裡包含了12個變量，其中特徵變量11個，結果變來呢為quality變量。該數據庫中將白酒品質總共分為1到10這10個等級，本文數據中包含了3至9這7個等級，為了方便分析，我們將白酒品質分為3個等級，其中品質3、4、5為“bad”品質，品質6為“mid”品質，品質7、8、9為“good”品質。分析代碼如下：

> wine <- read.csv("winequality-white.csv", sep = ";", header = T)

> summary(wine)

fixed.acidity volatile.acidity citric.acid residual.sugar chlorides free.sulfur.dioxide total.sulfur.dioxide

Min. : 3.800 Min. :0.0800 Min. :0.0000 Min. : 0.600 Min. :0.00900 Min. : 2.00 Min. : 9.0

1st Qu.: 6.300 1st Qu.:0.2100 1st Qu.:0.2700 1st Qu.: 1.700 1st Qu.:0.03600 1st Qu.: 23.00 1st Qu.:108.0

Median : 6.800 Median :0.2600 Median :0.3200 Median : 5.200 Median :0.04300 Median : 34.00 Median :134.0

Mean : 6.855 Mean :0.2782 Mean :0.3342 Mean : 6.391 Mean :0.04577 Mean : 35.31 Mean :138.4

3rd Qu.: 7.300 3rd Qu.:0.3200 3rd Qu.:0.3900 3rd Qu.: 9.900 3rd Qu.:0.05000 3rd Qu.: 46.00 3rd Qu.:167.0

Max. :14.200 Max. :1.1000 Max. :1.6600 Max. :65.800 Max. :0.34600 Max. :289.00 Max. :440.0

density pH sulphates alcohol quality

Min. :0.9871 Min. :2.720 Min. :0.2200 Min. : 8.00 Min. :3.000

1st Qu.:0.9917 1st Qu.:3.090 1st Qu.:0.4100 1st Qu.: 9.50 1st Qu.:5.000

Median :0.9937 Median :3.180 Median :0.4700 Median :10.40 Median :6.000

Mean :0.9940 Mean :3.188 Mean :0.4898 Mean :10.51 Mean :5.878

3rd Qu.:0.9961 3rd Qu.:3.280 3rd Qu.:0.5500 3rd Qu.:11.40 3rd Qu.:6.000

Max. :1.0390 Max. :3.820 Max. :1.0800 Max. :14.20 Max. :9.000

> #將白酒品質分為3個等級

> #設置中間變量對處理後的向量進行臨時存儲

> cha <- 0

> for(i in 1:4898) {

+ if (wine[i,12]>6) cha[i] <- "good"

+ else if(wine[i,12]>5) cha[i] <- "mid"

+ else cha[i] <- "bad"

+ }

> #將字符型變量轉化為含有因子的變量並賦值給數據集wine

> wine[,12]=factor(cha)

> summary(wine$quality)

bad good mid

1640 1060 2198

我們將利用數據集建立出適合的單隱藏層前饋人工神經網絡模型。在模型中我們將根據樣本白酒的非揮發性酸、揮發性酸、檸檬酸、剩餘糖分、氯化物、遊離二氧化硫、總二氧化硫、密度、酸性、硫酸鹽、酒精度這11個屬性來對白酒的品質進行判別。

數據處理

在建立人工神經網絡模型之前，我們首先應對數據進行預處理。

作為建立人工神經網絡模型的處理方式主要進行數據的歸一化。數據歸一化方法是神經網絡預測前對數據常做的一種處理方法，即將所有數據都轉化為[0,1]之間的束，其目的是取消各維度數據間數量級的差別，避免因為輸入輸出數據數量級別較大而造成網絡預測誤差較大。

數據歸一化的方法主要有以下兩種：

1、最大最小法。函數形式如下：

Xk = （Xk - Xmin）/（Xmax - Xmin）

其中Xmin為數據序列中的最小數，Xmax為序列中的最大數。

2、平均數方差法。函數形式如下：

Xk = （Xk - Xmean） / Xvar

其中，Xmean為數據序列的均值，Xvar為數據的方差。

下面採用最大最小法，對於這種0-1歸一化方法，我們將通過自寫程序對原始數據進行預處理，程序文件命名為“scale01.R”，相應代碼如下：

#確定程序名稱為scale01

scale01 <- function(x){

#提取預處理樣本集中特徵變量個數

ncol <- dim(x)[2]-1

#提取預處理樣本集中樣本總量

nrow <- dim(x)[1]

#建立用於保存新樣本集的矩陣

new <- matrix(0,nrow,ncol)

for(i in 1:ncol){

#提取每個變量的最大值

max <- max(x[,i])

#提取每個變量的最小值

min <- min(x[,i])

for(j in 1:nrow){

#計算歸一化後的新數據集

new[j,i] <- (x[j,i]-min)/(max-min)

}

}

new

}

建立模型

nnet()函數在建立支持單隱藏層前饋神經網絡模型的時候有兩種建立方式，一種是根據既定公式建立模型，而另一種是根據所給的數據建立模型。接下來我們將具體講述基於上述數據函數的兩種建模過程。

根據函數的第一種使用格式，在針對上述數據建模時，應該先確定我們所建立模型所使用的數據，然後再確定所建立模型的響應變量和自變量，具體建模操作如下：

> library(nnet)

> set.seed(71)

> #從總樣本集中抽取3000個樣本作為訓練集

> samp <- sample(1:4898,3000)

> #對樣本進行預處理

> source("scale01.R")

> wine[samp,1:11] <- scale(wine[samp,])

Error in colMeans(x, na.rm = TRUE) : 'x' must be numeric

> #確定參數rang的變化範圍

> r <- 1/max(abs(wine[samp,1:11]))

> set.seed(101)

> #建立神經網絡模型

> model1 <- nnet(quality~.,data=wine,subset=samp,size=4,rang=r,decay=5e-4,maxit=200)

在使用第一種格式建立模型時，如果使用數據中的全部自變量作為模型的自變量時，我們可以簡要的使用"quality~."代替全部的自變量。

根據函數的第二種使用格式，我們針對上述數據建立模型時，首先應該將響應變量和自變量分別提取出來。自變量通常用一個矩陣表示，而對於響應變量則應該進行相應的預處理。

首先要利用class.ind()函數將響應變量處理為類指標矩陣。在確定好數據後還應根據數據分析所使用的各項參數的具體值。對於建立神經網絡模型的具體過程如下：

> #提取wine數據集中除quality列以外的數據作為自變量

> x <- subset(wine,select = -quality)

> #提取wine數據集中的quality列數據作為響應變量

> y <- wine[,12]

> #對響應變量進行預處理

> y <- class.ind(y)

> set.seed(101)

> #建立神經網絡模型

> model2 <- nnet(x,y,decay = 5e-4,maxit = 200,size = 4,rang = r)

在使用第二種格式建立模型時，不需要特別強調所建立模型的形式，函數會自動將所有輸入到x矩陣中的數據作為建立模型所需要的自變量。在上述過程中，兩種模型的相關參數都是一樣的，兩種模型的權重衰減速度最小值為5e-4；最大迭代次數都為200次；隱藏層的節點數都為4個；最終我們建立出來的模型是一個11-4-3的神經網絡模型，即輸入層是11個節點，隱藏層是4個節點，輸出層是3個節點。

結果分析

我們使用summary()函數查看結果。

> summary(model1)

a 11-4-3 network with 63 weights

options were - softmax modelling decay=5e-04

b->h1 i1->h1 i2->h1 i3->h1 i4->h1 i5->h1 i6->h1 i7->h1 i8->h1 i9->h1 i10->h1 i11->h1

-52.03 -0.03 0.75 0.01 -0.02 0.39 0.00 0.00 53.67 -0.18 -0.24 -0.05

b->h2 i1->h2 i2->h2 i3->h2 i4->h2 i5->h2 i6->h2 i7->h2 i8->h2 i9->h2 i10->h2 i11->h2

2.42 -0.60 -4.12 2.12 -0.24 17.76 -0.06 0.02 7.41 -4.60 1.53 0.28

b->h3 i1->h3 i2->h3 i3->h3 i4->h3 i5->h3 i6->h3 i7->h3 i8->h3 i9->h3 i10->h3 i11->h3

-8.64 -11.41 0.41 -4.39 -17.36 0.07 -19.32 4.77 -8.63 0.38 13.47 -10.25

b->h4 i1->h4 i2->h4 i3->h4 i4->h4 i5->h4 i6->h4 i7->h4 i8->h4 i9->h4 i10->h4 i11->h4

22.81 16.27 -9.61 -66.33 -9.52 7.77 -5.45 0.51 22.76 9.51 -15.32 -6.63

b->o1 h1->o1 h2->o1 h3->o1 h4->o1

-13.58 27.08 6.70 -0.23 1.08

b->o2 h1->o2 h2->o2 h3->o2 h4->o2

12.61 -26.27 -16.02 0.75 -0.76

b->o3 h1->o3 h2->o3 h3->o3 h4->o3

0.97 -0.81 9.32 -0.53 -0.33

通過summary()函數我們可以得到關於模型的相關信息。在輸出結果的第一行我們可以看到模型的總體類型，該模型總共有三層，輸入層有11個節點，隱藏層有4個節點，輸出層有3個節點，該模型的權重總共有63個。

在輸出結果的第二行顯示的是模型中的相關參數的設定，在該模型的建立過程中，我們只設定了相應的模型權重衰減最小值，所以這裡顯示出了模型衰減最小值為5e-4。

接下來的第三部分是模型的具體判斷過程，其中的i系列代表的是輸入層的11個節點，h系列代表的是隱藏層的4個節點，而o系列代表的是輸出層的3個節點。對於b，我們可以將它理解為模型中的常數項。第三部分中的數字代表的是沒一個節點向下一個節點的輸入值的權重值。

預測判別

通常我們利用樣本數據建立模型之後，主要的目的都是利用模型來進行相應的預測和判別。在利用nnet()函數建立的模型進行預測時，我們將使用predict()函數對模型進行預測。

在使用predict()函數時，我們應該首先確認將要用於預測模型的類別。由於我們在建立模型時有兩種監理方式，而利用predict()函數進行預測的時候，對於兩種模型會存在兩種不同的預測結果，所以我們必須分清楚將要進行預測的模型是哪一類模型。具體操作如下：

針對第一種建模方式所建立的模型：

> #確認需要進行預測的樣本特徵矩陣

> x <- wine[,1:11]

> #根據模型model1對x數據進行預測

> pred <- predict(model1,x,type = "class")

> set.seed(110)

> #隨機挑選8個預測結果進行展示

> pred[sample(1:4898,8)]

[1] "bad" "mid" "good" "mid" "mid" "mid" "mid" "mid"

在進行數據預測時，我們主要注意的問題就是必須保證用於預測的自變量向量的個數同模型建立時使用的自變量向量個數一致，否則將無法預測結果。在使用predict()函數進行預測時，我們不用刻意去調整預測結果類型。通過上述預測結果的展示，我們可以看到predict()函數在預測時會自動識別預測結果的類型，並自動生成了相應的類別名稱。相對來說，利用第一種建模方式建立的模型在預測時較為方便。

針對第二種建模方式所建立的模型：

> #確認需要進行預測的樣本特徵矩陣

> xt <- wine[,1:11]

> #根據模型model2對xt數據進行預測

> pred <- predict(model2,xt)

> #查看預測結果維度

> dim(pred)

[1] 4898 3

> #隨機挑選4個預測結果進行展示

> pred[sample(1:4898,4),]

bad good mid

[1,] 0.7883393 0.010327155 0.2200430

[2,] 0.4667784 0.063874484 0.4121678

[3,] 0.1170516 0.415681237 0.5202564

[4,] 0.8431664 0.007020471 0.3159719

通過predict()函數對第二種模型進行預測，我們可以的看到預測結果是一個矩陣，而不像第一種模型那樣直接預測出模型中類別的名字。

在隨機挑選的4個預測結果中，我們可以看到每個樣本對應3種類別分別有3個數字，而這3個數字正是3個輸出結果的輸出值。這3個數的求和大約是等於1的，所以我們又可以將它簡要地看作概率，即樣本為其中某一類別的概率，對於樣本類別的判別則為概率最大的那一類。

因此對於上述預測結果我們需要將其進行進一步處理，處理之後才能直觀地看出樣本的預測類別。對於預測結果pred的處理，具體過程如下：

> #為3個類別確定名稱

> name <- c("bad","good","mid")

> #確定每行中最大值所在的列

> prednew <- max.col(pred)

> #根據預測結果將其變為相對應的類別名稱

> prednewn <- name[prednew]

> set.seed(201)

> #隨機挑選8個預測結果進行展示

> prednewn[sample(1:4898,8)]

[1] "mid" "bad" "mid" "bad" "mid" "mid" "mid" "bad"

通常在進行預測之後，我們還需要檢查模型預測的精度，這便需要用到table()函數對預測結果和真實結果做出對比展示。過程如下：

> #確定真實值的每行中最大值所在的列

> true <- max.col(y)

> #模型預測精度展示

> table(true,prednewn)

prednewn

true bad good mid

1 1088 8 544

2 80 175 805

3 535 93 1570

通過觀察table()函數對模型預測精度的展示結果，我們可以看到在模型預測時，模型將所有屬於屬於bad品質的白酒中的1088個樣本預測正確，但將另外8個樣本預測為good品質，並且將544個樣本預測為mid品質；模型將所有屬於good品質的白酒中的175個樣本預測正確，但將另外80個樣本預測為bad品質，並且將805個樣本預測為mid品質；模型將所有屬於mid品質的白酒中的1570個樣本預測正確，但將另外93個樣本預測為good品質，並且將535個樣本預測為bad品質。

模型差異分析

在利用nnet()函數建立模型的過程中，其中參數Wts的值我們通常默認為原始值。但是在nnet()函數中，參數Wts的值在建立模型的過程中用於迭代的權重初始值，該參數的默認值為系統隨機生成，換句話說，我們每次建立模型所使用的迭代初始值都是不相同的。因此我們再實際建模過程中會遇到這樣的現象：我們用同樣的數據，採取同樣的節點數，設定同樣的參數，但是最後會得到兩個不同的模型，甚至是兩個差異非常大的模型。

為了具體介紹該問題，我們先使用iris數據集進行舉例，首先我們利用下列語句建立模型model1以及模型model2，具體代碼如下：

>library(nnet)

>x <- subset(iris, select = -Species)

>y <- iris[,5]

>y <- class.ind(y)

>set.seed(101)

>#建立模型model1

>model1 <- nnet(x,y,range=1/max(abs(x)),size=4,maxit=500,decay=5e-4)

>#建立模型model2

>model2 <- nnet(x,y,range=1/max(abs(x)),size=4,maxit=500,decay=5e-4)

從建立模型的語句觀察，我們發現兩個模型應該是一樣的模型，但是通過對其進行具體分析，我們將發現兩個模型存在很大的差異，接下來我們從三個方面對模型差異進行分析。

1、模型是否因為迭代次數達到最大值而停止

如果模型的不同是因為建立模型時迭代次數達到最大值而停止迭代所導致的，那麼我們可以直接改變迭代的最大次數來使模型變得更加精確。具體查看方式如下：

>查看model1的迭代過程中是否達到迭代次數最大值

> model1$convergence

[1] 0

>查看model2的迭代過程中是否達到迭代次數最大值

> model2$convergence

[1] 0

從輸出結果中我們可以看到，兩個模型的迭代結果值都為0，這說明了再建立模型的過程中，迭代的停止並非是因為模型的迭代次數達到了最大迭代數。所以說明模型的最大迭代次數並不是影響兩個模型不同的主要原因。

2、模型迭代的最終值

模型迭代的最終值即為模型擬合標準同模型權重衰減值的和。在模型的輸出結果中，主要包含在模型的value中，該值越小說明模型擬合效果越好。我們對模型迭代的最終值的觀察過程及結果如下：

>#查看模型model1的迭代最終值

> model1$value

[1] 3.198636

>#查看模型model2的迭代最終值

> model2$value

[1] 2.598032

從輸出結果中我們可以看到，兩個模型的迭代最終值差異並不是很大。

因此對於因為初始迭代值不同而導致的模型不同的情況，我們可以使用該結果值來進行判斷，我們應該多運行幾次nnet()函數，而選擇所有模型中該結果值最小的一個模型作為最理想的模型。

3、觀察兩個模型的預測結果

人工神經網絡模型的預測效果是該模型最終最核心的作用，所以對於兩模型差異的情況，我們必須對模型的預測能力做出分析。

如果兩個模型在預測能力上顯示不出任何差異，那麼我們討論兩個模型不同也就失去了意義，因為我們所追求的就是模型的預測能力，所以在模型的差異問題上，我們最關心的也是兩個模型的預測能力的差異。觀察過程及結果如下：

> #為三個類別確定名稱

> name <- c("setosa","versicolor","virginica")

> #對模型model1進行預測

> pred1 <- name[max.col(predict(model1,x))]

> #對模型model2進行預測

> pred2 <- name[max.col(predict(model2,x))]

> table(iris$Species,pred1)

pred1

setosa versicolor virginica

setosa 50 0 0

versicolor 0 49 1

virginica 0 0 50

> table(iris$Species,pred2)

pred2

setosa versicolor virginica

setosa 50 0 0

versicolor 0 49 1

virginica 0 0 50

優化建模

在以上對nnet()函數的特別問題分析之後，我們瞭解到用相同數據相同參數建立的模型有可能不是最優模型。那麼，應該怎麼做才能得到最優模型呢？

針對這個問題，如果在時間和條件允許的情況下，我們可以多運行幾次模型，並從中挑選出針對於測試集樣本誤判率最小的模型。

首先，要確定出隱藏層最優節點的數目。之前已經介紹了對於人工神經網絡模型中隱藏層的相關確定條件，但是在實際模型構建過程中，仍需要儘可能地測試每一節點數目下模型的誤判率，以確定最優的模型誤判率。實現代碼如下：

> wine <- read.csv("winequality-white.csv", sep = ";", header = T)

> set.seed(71)

> wine <- wine[sample(1:4898,3000),]

> nrow.wine <- dim(wine)[1]

> source("scale01.R")

> #設置中間變量對處理後的向量進行臨時存儲

> cha <- 0

> for(i in 1:4898) {

+ if (wine[i,12]>6) cha[i] <- "good"

+ else if(wine[i,12]>5) cha[i] <- "mid"

+ else cha[i] <- "bad"

+ }

Error in if (wine[i, 12] > 6) cha[i] <- "good" else if (wine[i, 12] > :

missing value where TRUE/FALSE needed

> #將字符型變量轉化為含有因子的變量並賦值給數據集wine

> wine[,12]=factor(cha)

> set.seed(444)

> #從總樣本集中抽取70%的樣本作為訓練集

> samp <- sample(1:nrow.wine, nrow.wine*0.7)

> #對訓練集樣本進行預處理

> wine[samp,1:11] <- scale01(wine[samp,])

> wine[-samp,1:11] <- scale01(wine[-samp,])

> #確定參數range的變化範圍

> r <- 1/max(abs(wine[samp,1:11]))

> n <- length(samp)

> err1 <- 0

> err2 <- 0

> for(i in 1:17){

+ set.seed(111)

+ model <- nnet(quality~.,data=wine,maxit=400,rang=r,size=i,subset=samp,decay=5e-4)

+ err1[i] <- sum(predict(model,wine[samp,1:11],type="class")!=wine[samp,12])/n

+ err2[i] <- sum(predict(model,wine[-samp,1:11],type="class")!=wine[-samp,12])/(nrow.wine -n)

+ }

> plot(1:17,err1,'l',col=1,lty=1,ylab = "模型誤判率",xlab="隱藏層節點個數",ylim=c(min(min(err1),min(err2)),max(max(err1),max(err2))))

> lines(1:17,err2,col=1,lty=3)

> points(1:17,err1,col=1,pch="+")

> points(1:17,err2,col=1,pch="o")

> legend(1,0.53,"測試集誤判率",bty = "n",cex=1.5)

> legend(1,0.35,"訓練集誤判率",bty="n",cex=1.5)

經過上述程序運行之後，將得到關於樣本集在不同的隱藏層節點數下所對應的模型誤判率。從圖中我們可以清楚地看到，訓練集樣本錯誤跟隨隱藏層節點數的增加而下降，但是與此同時，測試集樣本錯誤卻未隨著隱藏層節點的增加而下降，這種現象便是由於模型中隱藏層節點數增加而引起的模型過度擬合導致的。

從圖中可以看到，模型針對測試集誤判率大概在模型隱藏層節點數為3的時候取到最小值，所以我們將隱藏層節點數數確定為3。

從前文中我們分析到，當神經網絡模型訓練週期過長的時候，建立的人工神經網絡模型將會記錄下訓練集中幾乎全部信息，這將會產生過度擬合的問題。即該模型針對於訓練集的時候將會體現出非常優異的預測能力，但是由於該模型記錄下了訓練集中的全部信息，則該模型也將訓練集中的許多特有的信息記錄下來，所以當模型用於其他樣本集的時候，模型的預測能力將會大大下降，即模型的泛化能力非常弱。

在確定最優隱藏層節點數的時候，接下來確定出最優的迭代次數，實現代碼如下：

>err11 <- 0

>err12 <- 0

>for(i in 1:500){

+set.seed(111)

+model <- nnet(quality~.,data=wine,maxit=i,rang=r,size=3,subset=samp)

+err11[i] <- sum(predict(model,wine[samp,1:11],type="class")!=wine[samp,12])/n

+err12[i] <- sum(predict(model,wine[-samp,1:11],type="class")!=wine[-samp,12])/(nrow.wine -n)

>}

> plot(1:length(err11),err11,'l',col=1,ylab = "模型誤判率",xlab="訓練週期",ylim=c(min(min(err11),min(err12)),max(max(err11),max(err12))))

> lines(1:length(err11),err12,col=1,lty=3)

> legend(250,0.47,"測試集誤判率",bty = "n",cex=1.2)

> legend(250,0.425,"訓練集誤判率",bty="n",cex=1.2)

從上圖可以看到，模型針對於訓練集和測試集的誤判率均同時隨訓練週期的增大而降低，之前也討論到當模型訓練週期過長時，模型應該會出現過度擬合的問題，即在訓練週期達到一定程度時，測試集誤差將會反向變化，訓練集誤差將會隨著模型訓練週期的增大而增大。

對於這個問題，是用R語言進行模型構建時會經常遇到，但這並非說明理論出現了錯誤。對該問題進行進一步分析可以得知出現該問題存在著兩個原因。

首先，在nnet程序包中，函數在構建模型時將會設定一個條件值以避免函數進入死循環。即在默認情況下，當函數計算值變化為零時模型將會停止運轉，所以很多時候模型將不會運行到過高的訓練週期。

其次，由於訓練集樣本同測試集樣本的相似度過高，所以訓練集中的特徵同樣為測試集中的特徵，所以即使在過度擬合的情況下，所構建的模型同樣能很好地適用於與訓練集相似度很高的數據集。

儘管會出現上圖中的問題，但是該圖像仍然具有一定的參考價值。從圖中可以發現，訓練集誤差隨著訓練週期的增大而不斷減小；但是對於測試集，當訓練週期達到一定程度後，模型的誤差率將會趨於平穩，模型的誤判率將不再下降。所圖中的情況，我們綜合分析決定將模型的訓練週期確定為300。

因此，最終得出的模型為隱藏層節點數為3，訓練週期為300，對於最新抽取的樣本集中，在隨機數生成器初始值為111情況下的人工神經網絡模型。

>set.seed(111)

>model <- nnet(quality~.,data = wine,maxit=300,rang=r,size=3,subset = samp)

> #根據需要進行預測的樣本特徵矩陣

> x <- wine[-samp,1:11]

> #根據模型model對x數據進行預測

> pred <- predict(model,x,type = "class")

> table(wine[-samp,12],pred)

pred

bad good mid

bad 237 5 51

good 20 59 109

mid 177 35 207

相關閱讀：《R語言數據挖掘實踐——用R語言實現神經網絡》