R語言判別分析

阿新 • • 發佈：2019-01-10

自己整理編寫的R語言常用資料分析模型的模板，原檔案為Rmd格式，直接複製貼上過來，作為個人學習筆記儲存和分享。部分參考薛毅的《統計建模與R軟體》和《R語言實戰》

本文中分三個方法介紹判別分析，Bayes判別，距離判別，Fisher判別。前兩種判別方法都要考慮兩個、或多個總體協方差（這裡是算方差，方差是協方差的一種）相等或不等的情況，由var.equal=的邏輯引數表示，預設是FALSE，表示認為兩總體協方差不等。用樣本的協方差可以估計總體的協方差。在Bayes方法中我們把相等和不等的兩個結果都列了出來，距離判別裡我們預設兩總體協方差不等。事實上，一般使用時，我們都以兩總體的協方差不等作為標準來進行後續計算。

協方差計算公式： S(xy) = N,i ∑ [(xi-x均)(yi-y均)]/(N-1) ; 這裡要計算的總體協方差就是方差，為 S = N,i ∑ [(xi-x均)*(xi-x均)^T]/(N-1)

1. Bayes判別

Bayes判別是假定對研究物件已有一定的認識，這種認識常用先驗概率來描述，當取得樣本後，就可以用樣本來修正已有的先驗概率分佈，得出後驗概率分佈，現通過後驗概率分佈進行各種統計推斷.

* 1.1 兩總體判別的Bayes判別程式*

在程式中，輸入變數TrnX1，TrnX2表示X1類，X2類訓練樣本，其輸入格式是資料框，或矩陣（樣本按行輸入）.rate= L(1|2) / L(2|1)*p2/p1 ，預設值為1. TstX是待測樣本，其輸入格式是資料框，或矩陣（樣本按行輸入），或向量（一個待測樣本）.如果不輸入TstX（預設值），則待測樣本為兩個訓練樣本之和，即計算訓練樣本的回代情況.輸入變數var.equal是邏輯變數，var.equal=TRUE表示認為兩總體的協方差陣是相同的；否則（預設值）是不同的. 函式的輸出是由”1”,”2”構成的一維矩陣，”1”代表待測樣本屬於X1類，”2”代表待測樣本屬於X2類

discriminiant.bayes<-function
   (TrnX1, TrnX2, rate=1, TstX = NULL, var.equal = FALSE){
   if (is.null(TstX) == TRUE) TstX<-rbind(TrnX1,TrnX2)
   if (is.vector(TstX) == TRUE)  TstX<-t(as.matrix(TstX))
   else if (is.matrix(TstX) != TRUE)
      TstX<-as.matrix(TstX)
   if (is.matrix(TrnX1 
) != TRUE) TrnX1<-as.matrix(TrnX1)
   if (is.matrix(TrnX2) != TRUE) TrnX2<-as.matrix(TrnX2)

   nx<-nrow(TstX)
   blong<-matrix(rep(0, nx), nrow=1, byrow=TRUE, 
         dimnames=list("blong", 1:nx))
   mu1<-colMeans(TrnX1); mu2<-colMeans(TrnX2) 
   if (var.equal == TRUE  || var.equal == T){
      S<-var(rbind(TrnX1,TrnX2)); beta<-2*log(rate)
      w<-mahalanobis(TstX, mu2, S)-mahalanobis(TstX, mu1, S)
   }
   else{
      S1<-var(TrnX1); S2<-var(TrnX2)
      beta<-2*log(rate)+log(det(S1)/det(S2))
      w<-mahalanobis(TstX, mu2, S2)-mahalanobis(TstX, mu1, S1)
   }

   for (i in 1:nx){
      if (w[i]>beta)
          blong[i]<-1
      else
          blong[i]<-2
   }
   blong
}

* 1.2 例. 用1.1中Bayes判別程式對冠心病和正常人進行判別分類*

為研究舒張期血壓與血漿膽固醇對冠心病的作用，調查了50~59歲的女冠心病人15例和正常人16例。他們的舒張壓（x1）與血漿膽固醇（x2）列在下表中。試用判別分析法建立判別冠心病與正常人的判別函式

#它們的先驗概率pi分別用15/31、16/31來估計
#載入資料
TrnX1<-matrix(
   c(9.86, 13.33, 14.66, 9.33, 12.80, 10.66, 10.66, 13.33, 13.33, 13.33, 12.00, 14.66, 13.33, 12.80, 13.33, 
     5.18, 3.73, 3.89, 7.10, 5.49, 4.09, 4.45, 3.63, 5.96, 5.70, 6.19, 4.01, 4.01, 3.63, 5.96),
   ncol=2)
TrnX2<-matrix(
   c(10.66, 12.53, 13.33, 9.33, 10.66, 10.66, 9.33, 10.66, 10.66, 10.66, 10.40, 9.33, 10.66, 10.66, 11.20, 9.33,
     2.07, 4.45, 3.06, 3.94, 4.45, 4.92, 3.68, 2.77, 3.21, 5.02, 3.94, 4.92, 2.69, 2.43, 3.42, 3.63),
   ncol=2)
#輸入待測樣本TrnX
TstX<-matrix(
   c(9.06,13.00,12.66,9.00,12.12,11.66,12.11,12.63,8.33,11.12,
     5.68,3.43,2.82,6.86,5.19,2.17,4.12,3.26,2.91,4.22 ),
   ncol=2)
#載入兩總體的貝葉斯判別函式    注 把貝葉斯判別函式存在了wd中
source("discriminiant.bayes.R")
#### 總體協方差陣相同
discriminiant.bayes(TrnX1, TrnX2,rate=16/15,var.equal=TRUE)    #rate=L(1|2)/L(2|1)* p2/p1
discriminiant.bayes(TrnX1, TrnX2,TstX,rate=16/15,var.equal=TRUE)
#### 總體協方差陣不同
discriminiant.bayes(TrnX1, TrnX2,rate=16/15)
discriminiant.bayes(TrnX1, TrnX2,TstX,rate=16/15)

在判別函式的鑑定下，有5個被錯判，3個樣本從冠心病組被錯判為正常組，兩個樣本從正常組錯判到冠心病組，分別為1,6,7,17,18

待測樣本分類結果為：

                 1 2 3 4 5 6 7 8 9 10  （樣本協方差陣相同）
           blong 2 1 2 1 1 2 1 2 2  2

                 1 2 3 4 5 6 7 8 9 10  （樣本協方差陣不同）
           blong 1 1 2 1 1 2 1 2 2  2

* 1.3 用R包裡的lda函式進行判別，得到判別函式*

#也可以csv格式匯入資料，先將以上兩類資料存為csv格式文件
coronary_disease <- read.csv("coronary_disease.csv")
#把分組變數變為定性變數
group <- factor(coronary_disease$group)
#隨機抽取20個一般樣本做訓練樣本
#train <- sample(1:31,20)
#顯示訓練樣本中各類的比例
#table(group[train])
#group作為分組變數，X1,X2作為判別變數，使用訓練樣本生成判別函式，先驗概率各為50%
library(MASS)
Z <- lda(group~.,data = coronary_disease ,prior=c(1,1)/2) #subset=train
Z

判別函式是Z = -0.6379195 * X1 + -0.8001452 * X2

* 1.4 取原資料80%做train，20%做test，檢視結果*

coronary_disease <- read.csv("coronary_disease.csv")
#讀取行數
N = length(coronary_disease$group)      
#ind=1的是0.7概率出現的行，ind=2是0.3概率出現的行
ind=sample(2,N,replace=TRUE,prob=c(0.8,0.2))
#生成訓練集(這裡訓練集和測試集隨機設定為原資料集的80%,20%)
coronary_train <- coronary_disease[ind==1,]
#生成測試集
coronary_test <- coronary_disease[ind==2,]
#固定這個28分組
set.seed(7)
#可以看到這裡的train中，i個1組，j個2組，共i+j個數據，則test裡有(31-i-j個)資料
i <- table(coronary_train$group)[[1]]
j <- table(coronary_train$group)[[2]]
i;j

* 1.4.1 用1.1中Bayes判別程式對測試集進行判別分類*

#輸入train的組別1、2，TrnXi；輸入test
TrnX1 <- coronary_train[1:i,2:3]
TrnX2 <- coronary_train[(i+1):(i+j),2:3]
TstX <- coronary_test[,2:3]
source("discriminiant.bayes.R")
#### 總體協方差陣相同
##它們的先驗概率pi分別用i/24、j/24來估計
discriminant_train1 <- discriminiant.bayes(TrnX1, TrnX2,rate=j/i,var.equal=TRUE) ; discriminant_train1  #rate=L(1|2)/L(2|1)* p2/p1 
discriminant_test1 <- discriminiant.bayes(TrnX1, TrnX2,TstX,rate=j/i,var.equal=TRUE) ; discriminant_test1
#### 總體協方差陣不同
discriminant_train2 <- discriminiant.bayes(TrnX1, TrnX2,rate=j/i) ; discriminant_train2
discriminant_test2 <- discriminiant.bayes(TrnX1, TrnX2,TstX,rate=j/i) ; discriminant_test2

* 1.4.2 匯出結果到csv檔案中*

#在原始訓練集，測試集中加入一列Is，表示是訓練集還是測試集 
coronary_train$Is <- c(rep("train",(i+j))) ;coronary_test$Is <- c(rep("test",(31-i-j)))
#在原始訓練集，測試集中分別加入總體協方差相同或不同時候的判定組別Discriminant
coronary_train$Belong_VarSame <- discriminant_train1[1:(i+j)]
coronary_test$Belong_VarSame <- discriminant_test1[1:(31-i-j)]
coronary_train$Belong_VarDiff <- discriminant_train2[1:(i+j)]
coronary_test$Belong_VarDiff <- discriminant_test2[1:(31-i-j)]
#合成一張表
coronary_result <- rbind(coronary_train,coronary_test)
write.csv(coronary_result,"coronary_result.csv")

* 1.4.3 總體協方差相同或不同時，train和test的正確率*

#以train裡的Var相同情況為例：
true_value_VarSame_train = coronary_result[coronary_result$Is=="train",]$Belong_VarSame
predict_value_VarSame_trian = coronary_result[coronary_result$Is=="train",]$group
#計算模型精確度
error_VarSame_train = predict_value_VarSame_trian-true_value_VarSame_train
accuracy_VarSame_train = (nrow(coronary_train)-sum(abs(error_VarSame_train)))/nrow(coronary_train) #精確度--判斷正確的數量佔總數的比例
accuracy_VarSame_train

2. 距離判別（二分類與多分類）

* 2.1 兩分類的距離判別程式*

在程式中，輸入變數TrnX1，TrnX2表示X1類，X2類訓練樣本，其輸入格式是資料框，或矩陣（樣本按行輸入）.TstX是待測樣本，其輸入格式是資料框，或矩陣（樣本按行輸入），或向量（一個待測樣本）.如果不輸入TstX（預設值），則待測樣本為兩個訓練樣本之和，即計算訓練樣本的回代情況.輸入變數var.equal是邏輯變數，var.equal=TRUE表示認為兩總體的協方差陣是相同的；否則（預設值）是不同的. 函式的輸出是由”1”,”2”構成的一維矩陣，”1”代表待測樣本屬於X1類，”2”代表待測樣本屬於X2類

discriminiant.distance<-function
   (TrnX1, TrnX2, TstX = NULL, var.equal = FALSE){
   if (is.null(TstX) == TRUE) TstX<-rbind(TrnX1,TrnX2)
   if (is.vector(TstX) == TRUE)  TstX<-t(as.matrix(TstX))
   else if (is.matrix(TstX) != TRUE)
      TstX<-as.matrix(TstX)
   if (is.matrix(TrnX1) != TRUE) TrnX1<-as.matrix(TrnX1)
   if (is.matrix(TrnX2) != TRUE) TrnX2<-as.matrix(TrnX2)

   nx<-nrow(TstX)
   blong<-matrix(rep(0, nx), nrow=1, byrow=TRUE, 
         dimnames=list("blong", 1:nx))
   mu1<-colMeans(TrnX1); mu2<-colMeans(TrnX2) 
   if (var.equal == TRUE  || var.equal == T){
      S<-var(rbind(TrnX1,TrnX2))
      w<-mahalanobis(TstX, mu2, S)-mahalanobis(TstX, mu1, S)
   }
   else{
      S1<-var(TrnX1); S2<-var(TrnX2)
      w<-mahalanobis(TstX, mu2, S2)-mahalanobis(TstX, mu1, S1)
   }
   for (i in 1:nx){
      if (w[i]>0)
          blong[i]<-1
      else
          blong[i]<-2
   }
   blong
}

* 2.2 多分類的距離判別程式*

distinguish.distance<-function
   (TrnX, TrnG, TstX = NULL, var.equal = FALSE){
   if ( is.factor(TrnG) == FALSE){
       mx<-nrow(TrnX); mg<-nrow(TrnG)
       TrnX<-rbind(TrnX, TrnG)
       TrnG<-factor(rep(1:2, c(mx, mg)))
   }
   if (is.null(TstX) == TRUE) TstX<-TrnX
   if (is.vector(TstX) == TRUE)  TstX<-t(as.matrix(TstX))
   else if (is.matrix(TstX) != TRUE)
      TstX<-as.matrix(TstX)
   if (is.matrix(TrnX) != TRUE) TrnX<-as.matrix(TrnX)

   nx<-nrow(TstX)
   blong<-matrix(rep(0, nx), nrow=1, dimnames=list("blong", 1:nx))
   g<-length(levels(TrnG))
   mu<-matrix(0, nrow=g, ncol=ncol(TrnX))
   for (i in 1:g)
      mu[i,]<-colMeans(TrnX[TrnG==i,]) 
   D<-matrix(0, nrow=g, ncol=nx)
   if (var.equal == TRUE  || var.equal == T){
      for (i in 1:g)
         D[i,]<- mahalanobis(TstX, mu[i,], var(TrnX),tol=2e-21)
   }
   else{
      for (i in 1:g)
         D[i,]<- mahalanobis(TstX, mu[i,], var(TrnX[TrnG==i,]),tol=2e-21)
   }
   for (j in 1:nx){
      dmin<-Inf
      for (i in 1:g)
          if (D[i,j]<dmin){
             dmin<-D[i,j]; blong[j]<-i
      }
   }
   blong
}

* 2.3 例2. 多分類的距離判別例項，判別冠心病和正常人*

例. 從心電圖的5個不同指標（x1~x5）對健康人（C=1）、血管硬化症患者（C=2）和冠心病患者（C=3）的資料進行判別分析，使用典型判別法。

#用和張佳玉sas文件裡一樣的資料輸入
X <- data.frame(
   x1 = c(8.11,9.36,9.85,2.55,6.01,9.64,4.11,8.9,7.71,7.51,8.06,6.8,8.68,5.67,8.1,3.71,5.37,5.22,4.71,4.71,3.36,8.27), 
   x2 = c(251.01,185.39,249.58,137.13,231.34,231.38,260.25,259.51,273.84,303.59,231.03,308.9,258.69,255.54,476.69,316.12,274.57,330.34,331.47,352.5,347.31,189.59),
   x3 = c(13.23,9.02,15.61,9.21,14.27,13.03,14.72,14.16,16.01,19.14,14.41,15.11,14.02,15.13,7.38,17.12,16.75,18.19,21.26,20.79,17.9,12.74),
   x4 = c(5.46,5.66,6.06,6.11,5.21,4.86,5.36,4.91,5.15,5.7,5.72,5.52,4.79,4.97,5.32,6.04,4.98,4.96,4.3,5.07,4.65,5.46),
   x5 = c(7.31,5.99,6.11,4.35,8.79,8.53,10.03,9.79,8.79,8.53,6.15,8.49,7.16,9.43,11.32,8.17,9.67,9.61,13.72,11,11.19,6.94)
)
G <- factor(c("1","1","1","1","1","1","1","1","1","1","1","2","2","2","2","2","2","3","3","3","3","3"),levels = c("1","2","3"))
source("distinguish.distance.R")
distinguish.distance(X,G)   #預設協方差不同

#用薛毅書中例子裡的資料輸入
X1 <- iris[,1:4]
G1 <- gl(3,50)
source("distinguish.distance.R")
distinguish.distance(X1,G1)    #預設協方差不同

這裡用薛毅的資料時沒有問題，換成自己的資料會一直報錯，後來找到了原因。

distinguish.distance.R程式中有兩行程式碼是mahalanobis(TstX, mu[i,], var(TrnX[TrnG==i,]))，因為我們的資料太小，出現如下錯誤：
（Error in solve.default(cov, …) :
system is computationally singular: reciprocal condition number = 2.60317e-20）
mahalanobis()函式會將這個很小的值認為是0，無法繼續後面的運算，所以報錯

解決辦法是，加一個引數,tol=2e-21，這個值設定成比error裡的值小就行。

3. Fisher判別

在程式中，輸入變數TrnX1，TrnX2表示X1類，X2類訓練樣本，其輸入格式是資料框，或矩陣（樣本按行輸入）.TstX是待測樣本，其輸入格式是資料框，或矩陣（樣本按行輸入），或向量（一個待測樣本）.如果不輸入TstX（預設值），則待測樣本為兩個訓練樣本之和，即計算訓練樣本的回代情況. 函式的輸出是由”1”,”2”構成的一維矩陣，”1”代表待測樣本屬於X1類，”2”代表待測樣本屬於X2類

discriminiant.fisher<-function(TrnX1, TrnX2, TstX = NULL){
   if (is.null(TstX) == TRUE)    TstX<-rbind(TrnX1,TrnX2)
   if (is.vector(TstX) == TRUE)  TstX<-t(as.matrix(TstX))
   else if (is.matrix(TstX) != TRUE)
      TstX<-as.matrix(TstX)
   if (is.matrix(TrnX1) != TRUE)  TrnX1<-as.matrix(TrnX1)
   if (is.matrix(TrnX2) != TRUE)  TrnX2<-as.matrix(TrnX2)

   nx<-nrow(TstX)
   blong<-matrix(rep(0, nx), nrow=1, byrow=TRUE, 
         dimnames=list("blong", 1:nx))
   n1<-nrow(TrnX1); n2<-nrow(TrnX2) 
   mu1<-colMeans(TrnX1); mu2<-colMeans(TrnX2) 
   S<-(n1-1)*var(TrnX1)+(n2-1)*var(TrnX2) 
   mu<-n1/(n1+n2)*mu1+n2/(n1+n2)*mu2
   w<-(TstX-rep(1,nx) %o% mu) %*% solve(S, mu2-mu1);
   for (i in 1:nx){
       if (w[i]<=0)
           blong[i]<-1
       else
           blong[i]<-2
   }
   blong
}

* 3.1 例3. Fisher判別分析法建立心肌梗塞病和正常人的判別函式*

為研究心肌梗塞病的危險因素，某研究者考察了兩組人群（心梗組與正常組，兩組間受試物件的年齡與性別構成接近）有關的指標10多項，現取其中血脂方面的6項指標：tc（總膽固醇）、tg（甘油三酯）、hdlc（高密度脂蛋白膽固醇）、ldlc（低密度脂蛋白膽固醇）、apoa（載脂蛋白AI）、apob（載脂蛋白B）。指標的測定結果如下表所示，每組各取30例，試作判別分析。

classX1 <- data.frame(
    x1 = c(245,236,238,233,240,235,204,200,297,177,200,195,166,144,233,143,228,264,178,240,180,161,236,168,174,215,268,213,285,193
), 
    x2 = c(157,275,354,250,149,166,365,95,240,97,172,211,217,111,107,91,223,186,131,127,211,91,106,106,141,168,185,387,154,123),
    x3 = c(38,40,38,31,35,40,38,43,38,49,43,47,33,28,42,24,34,41,49,33,27,39,36,36,28,38,28,22,39,42),
    x4 = c(168,125,126,150,170,164,90,100,207,108,116,106,86,46,156,108,136,183,98,174,106,88,104,104,103,134,203,141,210,121),
    x5 = c(1.1,1.22,0.9,1.02,1.26,1.3,1.33,1.24,1.14,1.49,1.25,1.22,1.1,0.71,0.95,0.67,1.05,1.22,1.18,0.78,0.85,0.94,0.87,0.87,0.81,0.88,0.75,0.8,1.17,1.12),
    x6 = c(1.01,1.12,1.06,0.98,1.13,1.15,0.95,0.98,1.51,1.02,1.03,0.94,0.74,0.65,0.77,0.65,0.84,0.92,1.27,0.9,0.69,0.52,0.58,0.73,0.73,0.87,0.97,0.78,1.37,1)
)
classX2 <- data.frame(
    x1 = c(174,106,173,178,198,180,134,204,168,219,189,180,177,172,166,210,166,223,136,156,201,134,195,262,194,165,183,200,171,222), 
    x2 = c(140,110,82,100,112,114,60,118,80,157,158,90,227,55,217,166,217,186,72,107,117,58,93,257,171,70,249,191,309,350),
    x3 = c(47,52,53,43,53,48,36,63,52,28,43,59,75,51,33,42,33,73,67,45,45,60,51,62,42,36,44,58,52,13),
    x4 = c(120,40,103,117,123,110,84,119,90,142,115,102,64,102,86,130,86,113,46,106,147,65,141,142,114,110,88,100,51,57),
    x5 = c(0.84,1.08,0.97,0.98,0.98,1.02,0.98,1.02,1.07,1.02,0.92,1.32,1.4,1.31,1.1,1.28,1.1,1.62,1.45,0.93,1.06,1.03,1.22,1.56,1.11,1.22,1.12,1.61,1.37,0.36),
    x6 = c(0.57,0.87,0.66,0.65,0.72,0.8,0.58,0.84,0.8,0.83,0.8,0.9,0.99,0.97,0.74,1.02,0.74,0.98,0.84,0.74,0.85,0.54,0.72,0.8,0.71,0.96,0.96,0.77,0.69,1.39)
)
source("discriminiant.fisher.R")
discriminiant.fisher(classX1,classX2)

將訓練樣本回代進行判別，有10個點判錯，分別是10、19、22、24、40、45、46、47、56、60號樣本.

#粗略匯出以上三個例子裡的結果
write.csv(discriminiant.bayes(TrnX1, TrnX2,TstX,rate=16/15,var.equal=TRUE),"disease1.csv")
write.csv(distinguish.distance(X1,G1),"disease2.csv")
write.csv(discriminiant.fisher(classX1,classX2),"disease3.csv")

R裡還可以用knn來進行判別，百度有很多博文有詳細介紹，這裡不作討論。

R語言判別分析

1. Bayes判別

2. 距離判別（二分類與多分類）

3. Fisher判別

R語言判別分析

R語言判別分析小結

R語言統計分析技術研究——嶺回歸技術的原理和應用

R語言關聯分析之啤酒和尿布

R語言︱情感分析—詞典型代碼實踐（最基礎）（一）

R語言生存分析

R語言迴歸分析函式說明

R語言bootstrap分析（boot）

R語言-錯誤分析-Error in .Call.graphics(C_palette2, .Call(C_palette2, NULL)) : invalid graphics state

分享《R語言資料分析與挖掘實戰(張良均等)》中文PDF+原始碼

R語言探索性分析及plyr資料轉換包

R語言生存分析

R語言因子分析

R語言統計分析篇

R語言——資料分析的一把利劍

R語言--關聯分析（問卷調查）

【R語言資料分析】豆瓣電影R語言爬蟲和資料分析

R語言學習-第九課-判別分析

應用統計學與R語言實現學習筆記（十一）——判別分析

R語言分類演算法之線性判別分析(Linear Discriminant Analysis)

R語言判別分析

1. Bayes判別

2. 距離判別（二分類與多分類）

3. Fisher判別

相關推薦