使用R進行描述性統計分析

title: “使用R進行描述性統計分析”

output:

html_document:

###單組計算描述性統計量

在對資料進行統計分析之前，一般我們需要知道目前資料的描述統計情況，如集中趨勢、離散趨勢、分布形狀。如果資料表中有很多變量，在excel中計算這些統計量的話，要利用公式一個一個進行計算，比較麻煩，在R中可以通過一些簡單的函數進行計算。

以車輛路試資料集mtcars為例，這裡我們隻需要三列資料，每加侖行駛英裡數mpg馬力hp和車重wt。

data(mtcars)
head(mtcars)
myvars <- c("mpg", "hp","wt")
head(mtcars[myvars])
           

####方法1：**summary()**函數：

summary()

為R自帶的函數，可以傳回最大值、最小值、四分位數、平均數6個值。

summary(mtcars[myvars])
           

####方法2：**apply()或者sapply()函數，

可以計算所選擇的任意描述統計量，對于sapply()**函數，其格式為sapply(x, FUN, options)，其中x為資料框，FUN為選擇的任意函數，例如mean(), sd(), var(), min(), max(), median(), length(), range(), quantile(), fivenum()。得出的結果偏度為+0.61（右偏态），峰度為-0.37，說明較正态分布稍平，

mystats <- function(x, na.omit=FALSE){
  if(na.omit)
    x <- x[!is.na(x)]  #如果x是非空值，則指派給x
  m <- mean(x)
  n <- length(x)
  s <- sd(x)
  skew <- sum((x-m)^3/s^3)/n
  kurt <- sum((x-m)^4/s^4)/n-3
  return(c(n=n, mean=m, stdev=s, skew=skew, kurt=kurt))
}
myvars <- c("mpg", "hp", "wt")
sapply(mtcars[myvars], mystats)
           

計算資料的峰度和偏度還可以用e1071包中的**skewness（）和kurtosis（）**函數來計算。

library(e1071)
mystats <- function(x, na.omit=FALSE){
  if(na.omit)
    x <- x[!is.na(x)]  
  skew <- skewness(x)
  kurt <- kurtosis(x)
  return(c(skew=skew, kurt=kurt))
}
myvars <- c("mpg", "hp", "wt")
sapply(mtcars[myvars], mystats)
           

####方法3：Hmisc包中的

describe()

函數

describe()

相比于方法1，可以傳回變量和觀測的數目，缺失值、唯一值數目，以及5個最大的值和5個最小的值。

library(Hmisc)
describe(mtcars[myvars])
           

####方法3：pastecs包中的

stat.desc()

函數

基本公式

stat.desc(x, basic=TRUE, desc=TRUE, norm=FALSE, p=0.95)

。

basic=TRUE

、

desc=TRUE

都是預設值，

norm=TRUE

（FALSE為預設值）傳回正态分布統計量，包括偏度、峰度以及它們的限制性水準和Shapiro-Wilk正态檢驗結果。該函數還計算了平均數95%的置信區間。

library(pastecs)
stat.desc(mtcars[myvars], norm=TRUE)
           

####方法4：psych包中的

describe()

函數

除了常用的統計量外，可以計算非缺失值的數量、截尾平均數、絕對中位數等。如果截尾平均數和平均數相差太大，說明存在極端值，此時截尾平均數可以更好地反映資料的集中趨勢。

注意：hmis包和psych包都含有

describe()

函數，但最後載入的程式包優先，hmis包

describe()

函數會被掩蔽，如果想要運作先前載入包的函數，可以使用代碼

hmisc::describe()

，告訴計算機執行哪個函數。

library(psych)
describe(mtcars[myvars])
           

###分組計算描述性統計量

如果一個變量有多個水準，想計算每個水準資料的描述統計情況，或者有多個被試組，計算每組被試的資料情況，這時可以使用分組計算的方法計算各組資料的描述統計量。

####方法1：**aggregate()**函數

list(am=mtcars$am),am指定了組的列标簽。**aggregate()**函數一次隻能傳回一種統計量

dstats <- function(x)sapply(x, mystats)
myvars <- c("mpg", "hp", "wt")
aggregate(mtcars[myvars], by=list(am=mtcars$am), mean)
aggregate(mtcars[myvars], by=list(am=mtcars$am), sd)
           

####方法2：**by()**函數

可以指定任意函數傳回多個指定的統計量，其基本公式為by(data, INDICES, FUN)。data為一個資料框；INDICES為一個因子或因子組成的清單，定義分組；FUN為任意函數。這裡的dstats調用了前面的mystats函數。

myvars <- c("mpg","hp","wt")
mystats <- function(x, na.omit=FALSE){
  if(na.omit)
    x <- x[!is.na(x)]  
  m <- mean(x)
  sd <- sd(x)
  skew <- skewness(x)
  kurt <- kurtosis(x)
  return(c(mean=m, SD=sd, skew=skew, kurt=kurt))
}
dstats <- function(x)sapply(x,mystats)
by(mtcars[myvars], mtcars$am, dstats)
           

####方法3：doBy包中的**summaryBy()**函數

可以指定任意函數傳回多個統計量。

library(doBy)
summaryBy(mpg + hp +wt~am, data=mtcars, FUN=mystats)
           

####方法4：**describeBy()**函數

傳回多個統計量，但不可以指定傳回哪些統計量

myvars <- c("mpg","hp","wt")
describeBy(mtcars[myvars], list(am=mtcars$am))