R語言多元統計分析初探

# 讀取多元統計分析數據到R
wine<-read.table("http://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data", sep=",")
# 繪制多元統計數據
# 矩陣散點圖
# 一種常見的方法是使用散點圖畫出多元統計數據，展現出所有變量兩兩之間的散點圖。
# 我們可以使用R中的“car”包里的“scatterplotMatrix()”函數來實現。
library(car)
scatterplotMatrix(wine[2:6])
# 組群標注數據點的散點圖
plot(wine$V4,wine$V5)
text(wine$V4,wine$V5,wine$V1,cex=0.7,pos=4,col="red")


# 輪廓圖？
# 輪廓圖？ 另一種非常有用的圖表類型便是”輪廓圖”，它通過繪制出每個變量在樣本中的值，展示出每個變量的變化。 
# 下文的“makeProfilePlot()”函數可以繪制出輪廓圖。這個函數需要“RColorBrewer”庫。
makeProfilePlot<-function(mylist,names){
  require(RColorBrewer)
  # find out how many variables we want to include
  numvariables<-length(mylist)
  # choose 'numvariables' random colours
  colours<-brewer.pal(numvariables,"Set1")
  # find out the minimum and maximum values of the variables:
  mymin<-1e+20
  mymax<-1e-20
  for(i in 1:numvariables){
    vectori<-mylist`i`
    mini<-min(vectori)
    maxi<-max(vectori)
    if(mini<mymin) {mymin<-mini}
    if(maxi>mymax) {mymax<-maxi}
  }
  
  # plot the variables
  for(i in 1:numvariables){
    vectori<-mylist`i`
    namei<-names[i]
    colouri<-colours[i]
    
    if(i == 1) {plot(vectori,col=colouri,type="l",ylim=c(mymin,mymax))}
    else  {points(vectori,col=colouri,type="l")}
    
    lastxval<-length(vectori)
    lastyval<-vectori[length(vectori)]
    text((lastxval-10),(lastyval),namei,col="black",cex=0.6)
  }
}
# 例如，為了畫出葡萄酒樣本中前五種化學物質的輪廓圖（他們存儲在“wine”變量的V2，V2，V4，V5，V6列），我們輸入：
library(RColorBrewer)
names<-c("V2","V3","V4","V5","V6")
mylist<-list(wine$V2,wine$V3,wine$V4,wine$V5,wine$V6)
makeProfilePlot(mylist,names)


# 計算多元統計數據的概要統計量
# 另一件事便是你可能會想計算你的多元統計數據集中每一個變量的概要統計量，像均值、標準偏差之類。
sapply(wine[,2:14],mean)
sapply(wine[,2:14],sd)
# 我們可以通過標準化來使數據看起來更有意義，以使我們能清楚的比較這些變量。我們需要便準化每一個變量以便使他們樣本方差為1，樣本均值為0.


# 每組的均值與方差
# 通常感興趣于從一個特定樣本群體去計算其均值和標準偏差，例如，計算每一個品種葡萄酒樣本。葡萄酒品種被存儲在“wine”變量的“V1”列中。
# 為了僅提取2號品種的數據，我們輸入：
cultivar2wine<-wine[wine$V1==2,]
sapply(cultivar2wine[2:14],mean)
sapply(cultivar2wine[2:14],sd)
# 你也可以通過相似的方法計算1號品種樣本，或者是3號品種樣本的13種化學物質濃度的均值和標準偏差： 
# 然而，為了方便起見，你也許想通過以下的“printMeanAndSdByGroup()”函數一次性輸出數據集中分組數據的均值和標準偏差：
printMeanAndSdByGroup<-function(variables,groupvariable){
  # find the names of the variables
  variablenames<-c(names(groupvariable),names(as.data.frame(variables)))
  # within each group, find the mean of each variable
  groupvariable<-groupvariable[,1] #ensures groupvariable is not a list
  means<-aggregate(as.matrix(variables)~groupvariable,FUN=mean)
  names(means)<-variablenames
  print(paste("Mean:"))
  print(means)
  # within each group, find the standard deviation of each variable:
  sds<-aggregate(as.matrix(variables)~groupvariable,FUN=sd)
  names(sds)<-variablenames
  print(paste("Standard deviations:"))
  print(sds)
  # within each group, find the number of samples:
  samplesizes<-aggregate(as.matrix(variables)~groupvariable,FUN=length)
  names(samplesizes)<-variablenames
  print(paste("Sample sizes:"))
  print(samplesizes)
}
printMeanAndSdByGroup(wine[2:14],wine[1])
# 函數”printMeanAndSdByGroup()”將輸出分組樣本的數字。在本例中，我們可以看到品種1有59個樣本，品種2有71個樣本，品種3有48個樣本。


## 變量的組間方差和組內方差
# 如果我們想計算特定變量的組內方差（例如，計算特定化學物質的濃度），我們可以使用下述的“calWithinGroupsVariance()”函數：
calcWithinGroupsVariance<-function(variable,groupvariable){
  # find out how many values the group variable can take
  groupvariable2<-as.factor(groupvariable`1`)
  levels<-levels(groupvariable2)
  numlevels<-length(levels)
  # get the mean and standard deviation for each group:
  numtotal<-0
  denomtotal<-0
  for(i in 1:numlevels){
    leveli<-levels[i]
    levelidata<-variable[groupvariable==leveli,]
    levelilength<-length(levelidata)
    # get the mean and standard deviation for group i:
    meani<-mean(levelidata)
    sdi<-sd(levelidata)
    numi<-(levelilength-1)*(sdi*sdi)
    denomi<-levelilength
    numtotal<-numtotal+numi
    denomtotal<-denomtotal+denomi
  }
  # calculate the within-groups variance
  Vw<-numtotal/(denomtotal-numlevels)
  return(Vw)
}
# 例如，計算V2變量（第一種化學物質的濃度）的組內方差，我們輸入：
calcWithinGroupsVariance(wine[2],wine[1]) # [1] 0.2620525
# 我們可以通過下述的“calcBetweenGroupsVariance()”函數來計算特定變量（如V2）的組間方差：
calcBetweenGroupsVariance <- function(variable,groupvariable) {
  # find out how many values the group variable can take 
  groupvariable2 <- as.factor(groupvariable`1`) 
  levels <- levels(groupvariable2)
  numlevels <- length(levels) 
  # calculate the overall grand mean: 
  grandmean <- mean(variable[,1]) 
  # get the mean and standard deviation for each group: 
  numtotal <- 0
  denomtotal <- 0 
  for (i in 1:numlevels) 
  {
    leveli <- levels[i] 
    levelidata <- variable[groupvariable==leveli,] 
    levelilength <- length(levelidata) 
    # get the mean and standard deviation for group i:
    meani <- mean(levelidata) 
    sdi <- sd(levelidata) 
    numi <- levelilength * ((meani - grandmean)^2)
    denomi <- levelilength 
    numtotal <- numtotal + numi 
    denomtotal <- denomtotal + denomi
  } 
  # calculate the between-groups variance 
  Vb <- numtotal / (numlevels - 1) 
  Vb <- Vb`1`
  return(Vb) 
}
# 可以像這樣使用它計算V2的組間方差：
calcBetweenGroupsVariance(wine[2],wine[1]) # [1] 35.39742
# 我們可以通過變量的組間方差除以組內方差計算“separation”。因此，這個通過V2計算的這個間隔是：
calcBetweenGroupsVariance(wine[2],wine[1])/calcWithinGroupsVariance(wine[2],wine[1])
# 如果我們想通過多元統計數據的所有變量計算出間隔，你可以使用下述的“calcSeparations()”:
calcSeparations<-function(variables,groupvariable){
  # find out how many variables we have
  variables<-as.data.frame(variables)
  numvariables<-length(variables)
  # find the variable names
  variablenames<-colnames(variables)
  # calculate the separation for each variable
  for(i in 1:numvariables){
    variablei<-variables[i]
    variablename<-variablenames[i]
    Vw<-calcWithinGroupsVariance(variablei,groupvariable)
    Vb<-calcBetweenGroupsVariance(variablei,groupvariable)
    sep<-Vb/Vw
    print(paste("variable",variablename,"Vw=",Vw,"Vb=",Vb,"separation=",sep))
  }
}
# 例如，計算每一個變量的13種化學物質濃度的間隔，我們輸入：
calcSeparations(wine[2:14],wine[1])
# 因此，個體變量在組內（葡萄酒品種）的最大間隔是V2（間隔為233.0）。
# 正如我們將在下面討論的，線性判別分析（LDA）的目的是尋找一個個體變量的線性組合將令組內（這里是品種）實現最大的間隔。
# 這里希望能夠通過任何個體變量（暫時是V8的233.9）得到一個更好的間隔替代這個最優間隔。