R的敘述性統計

許sir

本小節的目的為統整統計學在敘述性統計在R語言上的寫法

資料來源 https://rpubs.com/skydome20/R-Note10-Missing_Value

[基本資料設定]

  • 我們採用內建的iris資料

    1. 花萼長度(Sepal Length)
    2. 花萼寬度(Sepal Width)
    3. 花瓣長度(Petal Length)
    4. 花瓣寬度(Petal Width)
    5. 類別(Species)：可分為Setosa，Versicolor和Virginica三個品種。
In [2]:
data(iris)
dim(iris) # 查看列數與欄數
  1. 150
  2. 5

[PART 1]. 資料概述

  • 在收到一筆資料時, 我們常常先不去碰所謂的"統計學"
  • 而是先摸摸看資料, 看看它長的怎樣？是圓的還是方的？資料屬性是什麼？
  • 這樣做的目的, 是預先可以想像一下, 後續可能可以做些什麼？該用哪些統計方法來看資料？
In [10]:
#看看IRIS資料集的頭尾
head(iris, 10) #取前10筆資料
tail(iris,10) #取後10筆資料
Sepal.LengthSepal.WidthPetal.LengthPetal.WidthSpecies
5.1 3.5 1.4 0.2 setosa
4.9 3.0 1.4 0.2 setosa
4.7 3.2 1.3 0.2 setosa
4.6 3.1 1.5 0.2 setosa
5.0 3.6 1.4 0.2 setosa
5.4 3.9 1.7 0.4 setosa
4.6 3.4 1.4 0.3 setosa
5.0 3.4 1.5 0.2 setosa
4.4 2.9 1.4 0.2 setosa
4.9 3.1 1.5 0.1 setosa
Sepal.LengthSepal.WidthPetal.LengthPetal.WidthSpecies
1416.7 3.1 5.6 2.4 virginica
1426.9 3.1 5.1 2.3 virginica
1435.8 2.7 5.1 1.9 virginica
1446.8 3.2 5.9 2.3 virginica
1456.7 3.3 5.7 2.5 virginica
1466.7 3.0 5.2 2.3 virginica
1476.3 2.5 5.0 1.9 virginica
1486.5 3.0 5.2 2.0 virginica
1496.2 3.4 5.4 2.3 virginica
1505.9 3.0 5.1 1.8 virginica

[PART 2]. 基本趨勢量數

2-1.集中趨勢量數

In [18]:
mean(iris$Sepal.Length)     #「花萼長度」的平均值
median(iris$Sepal.Length)   #「花萼長度」的中位數
max(iris$Sepal.Length)      #「花萼長度」中的最大值
min(iris$Sepal.Length)      #「花萼長度」中的最小值
sum(iris$Sepal.Length)      #「花萼長度」加總
5.84333333333333
5.8
7.9
4.3
876.5
In [21]:
#類別資料算平均，得到NA不是夢
mean(iris$Species)

Warning message in mean.default(iris$Species):
“argument is not numeric or logical: returning NA”
<NA>

2-2. 分散趨勢量數

In [19]:
var(iris$Sepal.Length)      #「花萼長度」的變異數
sd(iris$Sepal.Length)       #「花萼長度」的標準差
range(iris$Sepal.Length)    #「花萼長度」最小值和最大值(全距)
0.685693512304251
0.828066127977863
  1. 4.3
  2. 7.9
In [20]:
quantile(iris$Sepal.Length, probs=0.25)  # 第一四分位數 
quantile(iris$Sepal.Length, probs=0.75)  # 第三四分位數
25%: 5.1
75%: 6.4

2-3. 一個一個key好累, 有沒有一個指令就可以搞定的

In [22]:
summary(iris)

  Sepal.Length    Sepal.Width     Petal.Length    Petal.Width   
 Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100  
 1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300  
 Median :5.800   Median :3.000   Median :4.350   Median :1.300  
 Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199  
 3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800  
 Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500  
       Species  
 setosa    :50  
 versicolor:50  
 virginica :50
In [25]:
#可是summary 指令有些東西沒有呈現, 還是要自己算
sd(iris$Sepal.Length)  #標準差
0.828066127977863
In [26]:
n=table(iris$Sepal.Length)
as.numeric(names(n))[which.max(n)] #眾數
5

[PART 3].遺漏值處理

  • 不管你是發問卷還是使用公司的資料集
  • 你都有可能碰到"很雷" 的dada-set
  • 這時你必需處理missing data

  • 在處理遺漏值時，大多數的人都會「直接移除資料」或是用「平均值來填補遺漏值」，但這樣的做法並不推薦：前者會讓資料減少，後者不會產生任何資訊。

  • 因此在遺漏值處理的手法上，最推崇的就是「k-Nearest Neighbours」或「mice套件」來填補遺漏值。其中，mice的全名為Multivariate Imputation via Chained Equations。

In [29]:
#我們先使用iris的資料集，讓資料中隨機產生遺漏值
install.packages("missForest")
require(missForest) # prodNA() function

also installing the dependency ‘itertools’

Updating HTML index of packages in '.Library'
Making 'packages.html' ... done
Loading required package: missForest
Loading required package: randomForest
randomForest 4.6-12
Type rfNews() to see new features/changes/bug fixes.
Loading required package: foreach
Loading required package: itertools
Loading required package: iterators
In [32]:
data <- prodNA(iris, noNA = 0.1) ## 在iris資料內，隨機產生10%的遺漏值
head(data)
Sepal.LengthSepal.WidthPetal.LengthPetal.WidthSpecies
5.1 3.5 1.4 0.2 setosa
4.9 3.0 1.4 0.2 setosa
4.7 3.2 1.3 0.2 setosa
4.6 3.1 1.5 NA setosa
5.0 3.6 1.4 0.2 setosa
5.4 3.9 1.7 0.4 NA

3-1. 直接移除有遺漏值的資料

In [33]:
# 當一筆資料是完整的，回傳TRUE；當一筆資料有遺漏值，回傳FALSE
complete.cases(data)
  1. TRUE
  2. TRUE
  3. TRUE
  4. FALSE
  5. TRUE
  6. FALSE
  7. TRUE
  8. TRUE
  9. TRUE
  10. FALSE
  11. FALSE
  12. FALSE
  13. TRUE
  14. TRUE
  15. TRUE
  16. FALSE
  17. FALSE
  18. FALSE
  19. FALSE
  20. TRUE
  21. TRUE
  22. TRUE
  23. FALSE
  24. TRUE
  25. TRUE
  26. TRUE
  27. FALSE
  28. TRUE
  29. TRUE
  30. TRUE
  31. FALSE
  32. TRUE
  33. TRUE
  34. FALSE
  35. FALSE
  36. TRUE
  37. TRUE
  38. TRUE
  39. TRUE
  40. TRUE
  41. FALSE
  42. TRUE
  43. FALSE
  44. TRUE
  45. FALSE
  46. TRUE
  47. FALSE
  48. TRUE
  49. TRUE
  50. TRUE
  51. TRUE
  52. TRUE
  53. TRUE
  54. TRUE
  55. TRUE
  56. FALSE
  57. TRUE
  58. TRUE
  59. TRUE
  60. TRUE
  61. FALSE
  62. TRUE
  63. FALSE
  64. FALSE
  65. TRUE
  66. FALSE
  67. TRUE
  68. FALSE
  69. TRUE
  70. FALSE
  71. TRUE
  72. TRUE
  73. TRUE
  74. TRUE
  75. FALSE
  76. TRUE
  77. TRUE
  78. TRUE
  79. TRUE
  80. FALSE
  81. TRUE
  82. TRUE
  83. FALSE
  84. TRUE
  85. FALSE
  86. FALSE
  87. TRUE
  88. FALSE
  89. TRUE
  90. FALSE
  91. FALSE
  92. FALSE
  93. TRUE
  94. TRUE
  95. TRUE
  96. TRUE
  97. TRUE
  98. TRUE
  99. TRUE
  100. FALSE
  101. TRUE
  102. FALSE
  103. FALSE
  104. FALSE
  105. FALSE
  106. TRUE
  107. TRUE
  108. TRUE
  109. TRUE
  110. FALSE
  111. TRUE
  112. TRUE
  113. FALSE
  114. TRUE
  115. FALSE
  116. TRUE
  117. TRUE
  118. FALSE
  119. FALSE
  120. FALSE
  121. TRUE
  122. FALSE
  123. FALSE
  124. FALSE
  125. TRUE
  126. FALSE
  127. FALSE
  128. FALSE
  129. FALSE
  130. TRUE
  131. TRUE
  132. FALSE
  133. TRUE
  134. FALSE
  135. TRUE
  136. FALSE
  137. TRUE
  138. TRUE
  139. TRUE
  140. FALSE
  141. TRUE
  142. FALSE
  143. TRUE
  144. FALSE
  145. TRUE
  146. FALSE
  147. FALSE
  148. FALSE
  149. TRUE
  150. TRUE
In [36]:
# 移除有遺漏值的資料
rm.data <- data[complete.cases(data), ]
dim(rm.data) # 查看列數與欄數
rm.data
  1. 89
  2. 5
Sepal.LengthSepal.WidthPetal.LengthPetal.WidthSpecies
15.1 3.5 1.4 0.2 setosa
24.9 3.0 1.4 0.2 setosa
34.7 3.2 1.3 0.2 setosa
55.0 3.6 1.4 0.2 setosa
74.6 3.4 1.4 0.3 setosa
85.0 3.4 1.5 0.2 setosa
94.4 2.9 1.4 0.2 setosa
134.8 3.0 1.4 0.1 setosa
144.3 3.0 1.1 0.1 setosa
155.8 4.0 1.2 0.2 setosa
205.1 3.8 1.5 0.3 setosa
215.4 3.4 1.7 0.2 setosa
225.1 3.7 1.5 0.4 setosa
245.1 3.3 1.7 0.5 setosa
254.8 3.4 1.9 0.2 setosa
265.0 3.0 1.6 0.2 setosa
285.2 3.5 1.5 0.2 setosa
295.2 3.4 1.4 0.2 setosa
304.7 3.2 1.6 0.2 setosa
325.4 3.4 1.5 0.4 setosa
335.2 4.1 1.5 0.1 setosa
365.0 3.2 1.2 0.2 setosa
375.5 3.5 1.3 0.2 setosa
384.9 3.6 1.4 0.1 setosa
394.4 3.0 1.3 0.2 setosa
405.1 3.4 1.5 0.2 setosa
424.5 2.3 1.3 0.3 setosa
445.0 3.5 1.6 0.6 setosa
464.8 3.0 1.4 0.3 setosa
484.6 3.2 1.4 0.2 setosa
⋮⋮⋮⋮⋮⋮
945.0 2.3 3.3 1.0 versicolor
955.6 2.7 4.2 1.3 versicolor
965.7 3.0 4.2 1.2 versicolor
975.7 2.9 4.2 1.3 versicolor
986.2 2.9 4.3 1.3 versicolor
995.1 2.5 3.0 1.1 versicolor
1016.3 3.3 6.0 2.5 virginica
1067.6 3.0 6.6 2.1 virginica
1074.9 2.5 4.5 1.7 virginica
1087.3 2.9 6.3 1.8 virginica
1096.7 2.5 5.8 1.8 virginica
1116.5 3.2 5.1 2.0 virginica
1126.4 2.7 5.3 1.9 virginica
1145.7 2.5 5.0 2.0 virginica
1166.4 3.2 5.3 2.3 virginica
1176.5 3.0 5.5 1.8 virginica
1216.9 3.2 5.7 2.3 virginica
1256.7 3.3 5.7 2.1 virginica
1307.2 3.0 5.8 1.6 virginica
1317.4 2.8 6.1 1.9 virginica
1336.4 2.8 5.6 2.2 virginica
1356.1 2.6 5.6 1.4 virginica
1376.3 3.4 5.6 2.4 virginica
1386.4 3.1 5.5 1.8 virginica
1396.0 3.0 4.8 1.8 virginica
1416.7 3.1 5.6 2.4 virginica
1435.8 2.7 5.1 1.9 virginica
1456.7 3.3 5.7 2.5 virginica
1496.2 3.4 5.4 2.3 virginica
1505.9 3.0 5.1 1.8 virginica

3-2. 用「平均數」來填補遺漏值

  • 直接刪除遺漏值不太好，因為會造成資訊損失(information loss)。
  • 所以我們常會採取「填補遺漏值」的手法
In [38]:
mean.data <- data

mean.1 <- mean(mean.data[, 1], na.rm = T)  # 取第一欄位的平均數
na.rows <- is.na(mean.data[, 1])           # 第一欄位中，有遺漏值存在的資料

mean.data[na.rows, 1] <- mean.1 # 用第一欄位的平均數，填補第一欄位的遺漏值
In [39]:
# 當一筆資料是完整的，回傳TRUE；當一筆資料有遺漏值，回傳FALSE
complete.cases(mean.data$Sepal.Length) #看一下還有沒有遺漏值
  1. TRUE
  2. TRUE
  3. TRUE
  4. TRUE
  5. TRUE
  6. TRUE
  7. TRUE
  8. TRUE
  9. TRUE
  10. TRUE
  11. TRUE
  12. TRUE
  13. TRUE
  14. TRUE
  15. TRUE
  16. TRUE
  17. TRUE
  18. TRUE
  19. TRUE
  20. TRUE
  21. TRUE
  22. TRUE
  23. TRUE
  24. TRUE
  25. TRUE
  26. TRUE
  27. TRUE
  28. TRUE
  29. TRUE
  30. TRUE
  31. TRUE
  32. TRUE
  33. TRUE
  34. TRUE
  35. TRUE
  36. TRUE
  37. TRUE
  38. TRUE
  39. TRUE
  40. TRUE
  41. TRUE
  42. TRUE
  43. TRUE
  44. TRUE
  45. TRUE
  46. TRUE
  47. TRUE
  48. TRUE
  49. TRUE
  50. TRUE
  51. TRUE
  52. TRUE
  53. TRUE
  54. TRUE
  55. TRUE
  56. TRUE
  57. TRUE
  58. TRUE
  59. TRUE
  60. TRUE
  61. TRUE
  62. TRUE
  63. TRUE
  64. TRUE
  65. TRUE
  66. TRUE
  67. TRUE
  68. TRUE
  69. TRUE
  70. TRUE
  71. TRUE
  72. TRUE
  73. TRUE
  74. TRUE
  75. TRUE
  76. TRUE
  77. TRUE
  78. TRUE
  79. TRUE
  80. TRUE
  81. TRUE
  82. TRUE
  83. TRUE
  84. TRUE
  85. TRUE
  86. TRUE
  87. TRUE
  88. TRUE
  89. TRUE
  90. TRUE
  91. TRUE
  92. TRUE
  93. TRUE
  94. TRUE
  95. TRUE
  96. TRUE
  97. TRUE
  98. TRUE
  99. TRUE
  100. TRUE
  101. TRUE
  102. TRUE
  103. TRUE
  104. TRUE
  105. TRUE
  106. TRUE
  107. TRUE
  108. TRUE
  109. TRUE
  110. TRUE
  111. TRUE
  112. TRUE
  113. TRUE
  114. TRUE
  115. TRUE
  116. TRUE
  117. TRUE
  118. TRUE
  119. TRUE
  120. TRUE
  121. TRUE
  122. TRUE
  123. TRUE
  124. TRUE
  125. TRUE
  126. TRUE
  127. TRUE
  128. TRUE
  129. TRUE
  130. TRUE
  131. TRUE
  132. TRUE
  133. TRUE
  134. TRUE
  135. TRUE
  136. TRUE
  137. TRUE
  138. TRUE
  139. TRUE
  140. TRUE
  141. TRUE
  142. TRUE
  143. TRUE
  144. TRUE
  145. TRUE
  146. TRUE
  147. TRUE
  148. TRUE
  149. TRUE
  150. TRUE

3-3.[進階方法 A]

[Tutorial on 5 Powerful R Packages used for imputing missing values] https://www.analyticsvidhya.com/blog/2016/03/tutorial-powerful-packages-imputing-missing-values/

3-4.[進階方法 B] 用K-Nearest Neighbours填補遺漏值

  • K-Nearest Neighbours(KNN)運用在遺漏值填補上的想法很簡單：

    1. 現在有一群學生的成績，包含國文、數學、自然，但老師不小心弄丟小明的國文考卷，於是小明的「國文」分數是遺漏值。
    2. 如果在不重考的狀況下，我們要給小明一個分數，該怎麼做？
    3. KNN的概念告訴我們，應該先看小明「數學和自然」的分數，看和哪些同學(K位)很相近，然後再拿那些同學(K位)的國文分數，取平均或加權平均(或是其他手法)後，當作小明的分數來填補。
    4. 一句話概括：「就是找和自己很像的K個鄰居，然後從他們身上複製自己所沒有的東西。」

這就是用KNN來填補遺漏值的想法。

install.packages("bitops"); install.packages("DMwR"); library(DMwR)

imputeData <- knnImputation(data); head(imputeData)

[PART 4]. 常見的統計圖表

In [8]:
h1<-hist(iris$Sepal.Length, breaks=10, xlab="Sepal Length", main="Histogram of Sepal Length") 
# break 切成 10組
# 但你會發現怎麼只取8組
# 別忘了最左最右的兩刀
In [9]:
h1$breaks #看一下它怎麼切的
  1. 4
  2. 4.5
  3. 5
  4. 5.5
  5. 6
  6. 6.5
  7. 7
  8. 7.5
  9. 8
In [10]:
seq(4,8,by=0.4) # 自己決定切割的組距 （從4到8 每0.4 切一組距）
  1. 4
  2. 4.4
  3. 4.8
  4. 5.2
  5. 5.6
  6. 6
  7. 6.4
  8. 6.8
  9. 7.2
  10. 7.6
  11. 8
In [11]:
h2<-hist(iris$Sepal.Length, breaks=seq(4,8,by=0.4))

4-2. 直方圖-常態分配

在直方圖中, 我們很常搭配常態分配曲線, 來看資料符不符合常態

In [12]:
h3<-hist(iris$Sepal.Length, breaks=seq(4,8,by=0.4), prob=TRUE) 
m = mean(iris$Sepal.Length)
std = sqrt(var(iris$Sepal.Length))

curve(dnorm(x, mean=m, sd=std), col="darkblue", lwd=2, add=TRUE) #add=TRUE 是將曲線疊加上去

4-4. 長條圖Bar Chart

In [13]:
plot(iris$Species, xlab="Species")

[補充說明] 直方圖 vs 長條圖

  • 都用來表示次數分配
    1. 長條圖的每一根柱子之間是分開的，而直方圖則是相連
    2. 長條圖的X軸用來描述離散變數，直方圖則描述連續變數
      • nominal variable或ordinal variable適用長條圖
      • Interval variable或ratio variable則是用直方圖

4-5. 盒鬚圖Box Plot

In [25]:
par(mfrow=c(1,2))
boxplot(iris[,1:2]) #單變數畫圖
boxplot(iris$Sepal.Length~iris$Species) #兩個變數

dev.off()
null device: 1
In [22]:
par(mfrow=c(2,2))
boxplot(Sepal.Length  ~ Species, iris, main = "Sepal Length wrt Species", col = "lightpink3")
boxplot(Sepal.Width   ~ Species, iris, main = "Sepal Width wrt Species", col = "antiquewhite1")
boxplot(Petal.Length  ~ Species, iris, main = "Petal Length wrt Species", col = "lightskyblue4")
boxplot(Petal.Width  ~ Species, iris, main = "Petal Width wrt Species", col = "orange1")

dev.off()

4-6. 散布圖 Scatter Plot

In [16]:
plot(iris$Petal.Length, iris$Petal.Width, pch=19, col=iris$Species)
abline(lm(iris$Petal.Width~iris$Petal.Length), col="blue") #abline是畫一條線（本例而言是畫迴歸線 lm)

[散布矩陣圖] 用意

  • 觀察變數間的相關性
  • 變數的篩選：

  • 相關係數(Correlation Coefficient)

    #### [補充]：

    • 主成分分析(Principal Component Analysis)
    • Step-wise forward selection
    • Step-wise backward selection
In [19]:
with(iris, plot(iris, col = Species)) #散布矩陣圖

#你會看到最右邊與最下面呈現是名目資料
# 另外 Petal.Length 及 Petal.Width 呈現線性關系 （在我還沒有用其它任何統計方法, 我就大概猜的到這兩個變數有關係）

4-7. 圓餅圖 Pie chat

In [20]:
table(iris$Species)
pie(table(iris$Species), main = "Pie Chart of the Iris data set Species", 
    col = c("orange1", "chocolate", "coral"), radius = 1)

    setosa versicolor  virginica 
        50         50         50
In [ ]: