R語言多元統(tǒng)計(jì)分析初探

# 讀取多元統(tǒng)計(jì)分析數(shù)據(jù)到R
wine<-read.table("http://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data", sep=",")
# 繪制多元統(tǒng)計(jì)數(shù)據(jù)
# 矩陣散點(diǎn)圖
# 一種常見的方法是使用散點(diǎn)圖畫出多元統(tǒng)計(jì)數(shù)據(jù)，展現(xiàn)出所有變量?jī)蓛芍g的散點(diǎn)圖。
# 我們可以使用R中的“car”包里的“scatterplotMatrix()”函數(shù)來實(shí)現(xiàn)。
library(car)
scatterplotMatrix(wine[2:6])
# 組群標(biāo)注數(shù)據(jù)點(diǎn)的散點(diǎn)圖
plot(wine$V4,wine$V5)
text(wine$V4,wine$V5,wine$V1,cex=0.7,pos=4,col="red")


# 輪廓圖？
# 輪廓圖？ 另一種非常有用的圖表類型便是”輪廓圖”，它通過繪制出每個(gè)變量在樣本中的值，展示出每個(gè)變量的變化。 
# 下文的“makeProfilePlot()”函數(shù)可以繪制出輪廓圖。這個(gè)函數(shù)需要“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)
  }
}
# 例如，為了畫出葡萄酒樣本中前五種化學(xué)物質(zhì)的輪廓圖（他們存儲(chǔ)在“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)


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


# 每組的均值與方差
# 通常感興趣于從一個(gè)特定樣本群體去計(jì)算其均值和標(biāo)準(zhǔn)偏差，例如，計(jì)算每一個(gè)品種葡萄酒樣本。葡萄酒品種被存儲(chǔ)在“wine”變量的“V1”列中。
# 為了僅提取2號(hào)品種的數(shù)據(jù)，我們輸入：
cultivar2wine<-wine[wine$V1==2,]
sapply(cultivar2wine[2:14],mean)
sapply(cultivar2wine[2:14],sd)
# 你也可以通過相似的方法計(jì)算1號(hào)品種樣本，或者是3號(hào)品種樣本的13種化學(xué)物質(zhì)濃度的均值和標(biāo)準(zhǔn)偏差： 
# 然而，為了方便起見，你也許想通過以下的“printMeanAndSdByGroup()”函數(shù)一次性輸出數(shù)據(jù)集中分組數(shù)據(jù)的均值和標(biāo)準(zhǔn)偏差：
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])
# 函數(shù)”printMeanAndSdByGroup()”將輸出分組樣本的數(shù)字。在本例中，我們可以看到品種1有59個(gè)樣本，品種2有71個(gè)樣本，品種3有48個(gè)樣本。


## 變量的組間方差和組內(nèi)方差
# 如果我們想計(jì)算特定變量的組內(nèi)方差（例如，計(jì)算特定化學(xué)物質(zhì)的濃度），我們可以使用下述的“calWithinGroupsVariance()”函數(shù)：
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)
}
# 例如，計(jì)算V2變量（第一種化學(xué)物質(zhì)的濃度）的組內(nèi)方差，我們輸入：
calcWithinGroupsVariance(wine[2],wine[1]) # [1] 0.2620525
# 我們可以通過下述的“calcBetweenGroupsVariance()”函數(shù)來計(jì)算特定變量（如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) 
}
# 可以像這樣使用它計(jì)算V2的組間方差：
calcBetweenGroupsVariance(wine[2],wine[1]) # [1] 35.39742
# 我們可以通過變量的組間方差除以組內(nèi)方差計(jì)算“separation”。因此，這個(gè)通過V2計(jì)算的這個(gè)間隔是：
calcBetweenGroupsVariance(wine[2],wine[1])/calcWithinGroupsVariance(wine[2],wine[1])
# 如果我們想通過多元統(tǒng)計(jì)數(shù)據(jù)的所有變量計(jì)算出間隔，你可以使用下述的“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))
  }
}
# 例如，計(jì)算每一個(gè)變量的13種化學(xué)物質(zhì)濃度的間隔，我們輸入：
calcSeparations(wine[2:14],wine[1])
# 因此，個(gè)體變量在組內(nèi)（葡萄酒品種）的最大間隔是V2（間隔為233.0）。
# 正如我們將在下面討論的，線性判別分析（LDA）的目的是尋找一個(gè)個(gè)體變量的線性組合將令組內(nèi)（這里是品種）實(shí)現(xiàn)最大的間隔。
# 這里希望能夠通過任何個(gè)體變量（暫時(shí)是V8的233.9）得到一個(gè)更好的間隔替代這個(gè)最優(yōu)間隔。