機器學習及演算法
林嶔 (Lin, Chin)
Lesson 3 統計分析實作2（單變項統計分析）
本節課程目標介紹
– 請至這裡下載範例資料
dat <- read.csv("ECG_train.csv", header = TRUE, fileEncoding = 'CP950', stringsAsFactors = FALSE, na.strings = "")
- 今天的重點在於讓你進一步加深對R語言的依賴程度，我們期望用他來解決你日常的分析工作。
– 我們今天要教的檢定方式如下，這大約可以涵蓋90%以上的常見狀況：
第一節：列聯表統計(1)
- 我們先從最簡單的列聯表統計開始，若兩變項皆為類別變項，我們最常使用卡方檢定及Fisher exact test來進行檢定
– 這時候我們需要使用函數「chisq.test()」及「fisher.test()」
- 需要注意的是，上述這兩個方法都是對『列聯表』做計算，因此我們需要先使用函數「table()」幫助我們產生列聯表
TABLE <- table(dat[,"GENDER"], dat[,"AMI"])
TABLE
##
## not-AMI NSTEMI STEMI
## female 521 32 36
## male 539 136 219
Result1 <- chisq.test(TABLE)
Result1
##
## Pearson's Chi-squared test
##
## data: TABLE
## X-squared = 139.18, df = 2, p-value < 2.2e-16
Result2 <- fisher.test(TABLE)
Result2
##
## Fisher's Exact Test for Count Data
##
## data: TABLE
## p-value < 2.2e-16
## alternative hypothesis: two.sided
第一節：列聯表統計(2)
- 還記得上節課教的嗎?物件『Result1』以及物件『Result2』都是基於列表的物件，因此我如果最關心的是p-value，我們可以透過函數「ls()」看看這兩個物件裡面有甚麼，是否有p-value
## [1] "data.name" "expected" "method" "observed" "parameter" "p.value"
## [7] "residuals" "statistic" "stdres"
## [1] "alternative" "data.name" "method" "p.value"
- 透過下面的方法你就可以把p-value提取出來了：
## [1] 6.004532e-31
## [1] 1.192588e-33
第一節：列聯表統計(3)
- 判斷何時該使用卡方檢定或是Fisher exact test是比較困難的，目前的可以接受的說法是『有超過80%的格子期望值大於5』就使用卡方檢定，否則使用Fisher exact test
– 如果你有仔細注意的話，你會發現物件『Result1』裡面似乎有一個叫做「expected」的東西，而這個其實就是『有超過80%的格子期望值大於5』中，所謂的期望值：
##
## not-AMI NSTEMI STEMI
## female 420.998 66.72421 101.2778
## male 639.002 101.27579 153.7222
##
## not-AMI NSTEMI STEMI
## female TRUE TRUE TRUE
## male TRUE TRUE TRUE
mean(Result1$expected > 5)
## [1] 1
第一節：列聯表統計(4)
- 我們來試試看在一個子資料集中計算看看好了，這樣樣本比較容易不夠：
subdat <- dat[dat$rhythm.3 %in% 1,]
TABLE <- table(subdat[,"GENDER"], subdat[,"AMI"])
TABLE
##
## not-AMI NSTEMI STEMI
## female 6 0 0
## male 7 2 3
Result1 <- chisq.test(TABLE)
Result1$expected
##
## not-AMI NSTEMI STEMI
## female 4.333333 0.6666667 1
## male 8.666667 1.3333333 2
mean(Result1$expected > 5)
## [1] 0.1666667
- 看起來是沒有滿足條件，所以這個狀況下我們需要執行Fisher exact test
Result2 <- fisher.test(TABLE)
Result2
##
## Fisher's Exact Test for Count Data
##
## data: TABLE
## p-value = 0.3067
## alternative hypothesis: two.sided
練習1：計算卡方檢定的期望值
##
## not-AMI NSTEMI STEMI
## female 4.333333 0.6666667 1
## male 8.666667 1.3333333 2
##
## not-AMI NSTEMI STEMI
## female 6 0 0
## male 7 2 3
練習1答案(1)
- 我們首先需要計算Row and column marginals:
row_m <- c(sum(TABLE[1,]), sum(TABLE[2,]))
row_m
## [1] 6 12
col_m <- c(sum(TABLE[,1]), sum(TABLE[,2]), sum(TABLE[,3]))
col_m
## [1] 13 2 3
## [1] 18
- 有了這些工具後，我們就能進行運算了，我們先計算row=1&column=1的，依此類推：
row_m[1] * col_m[1] / size
## [1] 4.333333
練習1答案(2)
- 這個過程其實還可以再快一點，我們現在介紹函數「apply()」，他可以用來進行快速大量的運算：
row_m <- apply(TABLE, MARGIN = 1, FUN = sum)
row_m
## female male
## 6 12
col_m <- apply(TABLE, MARGIN = 2, FUN = sum)
col_m
## not-AMI NSTEMI STEMI
## 13 2 3
- 如果你會矩陣乘法，在這李我們介紹一下他的運算子「%*%」及函數「t()」進行轉置，你可以迅速的算出這個結果：
## not-AMI NSTEMI STEMI
## [1,] 13 2 3
row_m %*% t(col_m) / size
## not-AMI NSTEMI STEMI
## [1,] 4.333333 0.6666667 1
## [2,] 8.666667 1.3333333 2
第二節：獨立t檢定與M-W U檢定(1)
- 在自變項為2類別，依變相為連續變項時，我們最常使用獨立t檢定及M-W U檢定來進行檢定
– 其中M-W U檢定有另一個名字，叫做Wilcoxon rank sum test
– 因此他們相對應的函數是「t.test()」及「wilcox.test()」
- 我們使用性別做分組，檢定他們鉀離子(K)的差異，分別使用t test及Wilcoxon rank sum test
result1 <- t.test(dat[,"K"]~dat[,"GENDER"], var.equal = FALSE)
result1
##
## Welch Two Sample t-test
##
## data: dat[, "K"] by dat[, "GENDER"]
## t = -4.0971, df = 3147.5, p-value = 4.289e-05
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.28235100 -0.09957477
## sample estimates:
## mean in group female mean in group male
## 3.680688 3.871651
result2 <- wilcox.test(dat[,"K"]~dat[,"GENDER"])
result2
##
## Wilcoxon rank sum test with continuity correction
##
## data: dat[, "K"] by dat[, "GENDER"]
## W = 1108000, p-value = 1.459e-09
## alternative hypothesis: true location shift is not equal to 0
第二節：獨立t檢定與M-W U檢定(2)
- 知道為什麼直接列印出result1及result2會出現這麼多文字吧?因為函數「t.test()」及「wilcox.test()」輸出的結果跟之前遇到的差不多
## [1] "htest"
## [1] "htest"
- 當然，如果我們可以使用函數「ls()」看看這兩個物件裡面有多少東西
## [1] "alternative" "conf.int" "data.name" "estimate" "method"
## [6] "null.value" "parameter" "p.value" "statistic"
## [1] "alternative" "data.name" "method" "null.value" "parameter"
## [6] "p.value" "statistic"
- 如果我們最需要的是裡面的p value，各位知道要怎麼取得吧?
第二節：獨立t檢定與M-W U檢定(3)
- t test有分成兩種：『變異數同質』及『變異數不同質』。一般來說我們會先做變異數同值檢定，我們可以透過函數「var.test()」來做這個檢定
var.result <- var.test(dat[,"K"]~dat[,"GENDER"])
var.result
##
## F test to compare two variances
##
## data: dat[, "K"] by dat[, "GENDER"]
## F = 1.0081, num df = 1511, denom df = 1671, p-value = 0.8718
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.9136993 1.1125070
## sample estimates:
## ratio of variances
## 1.00809
- 看起來變異數是同質的，因此其實我們應該使用的是變異數同質的t test：
result3 <- t.test(dat[,"K"]~dat[,"GENDER"], var.equal = TRUE)
result3
##
## Two Sample t-test
##
## data: dat[, "K"] by dat[, "GENDER"]
## t = -4.0979, df = 3182, p-value = 4.273e-05
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.28233211 -0.09959366
## sample estimates:
## mean in group female mean in group male
## 3.680688 3.871651
第三節：ANOVA與K-W檢定(1)
- 若自變項為3類別以上時，依變相為連續變項時，我們通常會使用ANOVA and Kruskal-Wallis test來進行檢定
– 我們使用AMI做分組，檢定他們AGE的差異，分別使用ANOVA and Kruskal-Wallis test
- 我們需要使用函數「aov()」及「anova()」執行ANOVA
Variance.table <- aov(dat[,"AGE"]~as.factor(dat[,"AMI"]))
ANOVA.table <- anova(Variance.table)
ANOVA.table
## Analysis of Variance Table
##
## Response: dat[, "AGE"]
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(dat[, "AMI"]) 2 1296 647.93 2.0921 0.1238
## Residuals 1480 458360 309.70
- 函數「kruskal.test()」可以幫助我們做Kruskal-Wallis test
KW.result <- kruskal.test(dat[,"AGE"]~as.factor(dat[,"AMI"]))
KW.result
##
## Kruskal-Wallis rank sum test
##
## data: dat[, "AGE"] by as.factor(dat[, "AMI"])
## Kruskal-Wallis chi-squared = 8.7746, df = 2, p-value = 0.01243
- 聰明的你一定想的到，如果我們最關心的是p value，要怎樣拿到p value呢?
練習2：如何分辨何時要使用哪個檢定呢
- 對於「X = 類別變項」同時「Y = 連續變項」的情境，總共有4種檢定方式，我們希望你有能力區分何時該用什麼檢定
– 對於X而言，你必須有能力判斷他的類別數目
– 對於Y而言，你必須有能力判斷他在各組別的樣本數，如果有任一組樣本數<25那就必須使用相對應的無母數統計
練習2答案
- 對於第一個問題來說，你最需要想到的是利用第一節課教的因子向量進行操作，我們來看一下：
length(levels(as.factor(dat$GENDER)))
## [1] 2
length(levels(as.factor(dat$AMI)))
## [1] 3
- 計算樣本數乍看之下好像也不難，用函數「table()」好像可以快速的統計：
##
## female male
## 2172 2828
##
## not-AMI NSTEMI STEMI
## 1060 168 255
- 但這樣其實是有問題的，假設我們想要檢定的Y是K，他有部分NA必須排除，所以整個流程應該是這樣：
subdat <- dat[!dat$K %in% NA,]
table(subdat$GENDER)
##
## female male
## 1512 1672
##
## not-AMI NSTEMI STEMI
## 908 64 155
練習3：ANOVA的事後檢定
- 使用ANOVA檢定AMI與AGE的關係時，我們只能知道整體是否有差異，你是否有辦法做出下列三個比較：
「AMI = ‘STEMI’」比上「AMI = ‘NSTEMI’」
「AMI = ‘STEMI’」比上「AMI = ‘not-AMI’」
「AMI = ‘NSTEMI’」比上「AMI = ‘not-AMI’」
練習3答案
- 對於第一個問題來說，你需要將樣本進行限制，這樣就能執行「AMI = ‘STEMI’」比上「AMI = ‘NSTEMI’」的檢定了：
subdat <- dat[dat$AMI %in% c('STEMI', 'NSTEMI'),]
var.result <- var.test(subdat[,"AGE"]~subdat[,"AMI"])
var.result
##
## F test to compare two variances
##
## data: subdat[, "AGE"] by subdat[, "AMI"]
## F = 1.0207, num df = 167, denom df = 254, p-value = 0.8768
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.7768977 1.3518136
## sample estimates:
## ratio of variances
## 1.020698
t.result <- t.test(subdat[,"AGE"]~subdat[,"AMI"], var.equal = TRUE)
t.result
##
## Two Sample t-test
##
## data: subdat[, "AGE"] by subdat[, "AMI"]
## t = 1.0627, df = 421, p-value = 0.2885
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -1.215674 4.077034
## sample estimates:
## mean in group NSTEMI mean in group STEMI
## 62.61472 61.18404
– 他的邏輯是，如果檢定n次，那顯著水準就應該從本來的0.05改成0.05/n
– 換言之就是假定要保持顯著水準固定為0.05，那我們就必須將t檢定的結果乘上n
– 這是最終的p-value：
## [1] 0.8656402
第四節：相關性檢定(1)
- Pearson correlation是主要拿來做連續變項間相關性的統計方法，而Spearman correlation則是在樣本數不足時的替代方案
– 函數「cor.test()」可以用來計算相關係數
Result3 <- cor.test(dat[,"PR"], dat[,"K"], method = "pearson")
Result3
##
## Pearson's product-moment correlation
##
## data: dat[, "PR"] and dat[, "K"]
## t = 1.3747, df = 2820, p-value = 0.1693
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.01102990 0.06271683
## sample estimates:
## cor
## 0.02587868
Result4 <- cor.test(dat[,"PR"], dat[,"K"], method = "spearman")
Result4
##
## Spearman's rank correlation rho
##
## data: dat[, "PR"] and dat[, "K"]
## S = 3734100000, p-value = 0.8706
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.003068111
- 如何計算樣本是否足夠呢?你首先需要保證PR跟K都沒有NA：
sum(!(dat[,"PR"] %in% NA) & !(dat[,"K"] %in% NA))
## [1] 2822
第五節：相關係數(1)
- 在剛剛的檢定中你應該有發現，相關性檢定除了提供p-value外，更重要的是他有提供一個「相關係數」，這可以幫助我們量化關聯性的大小
– 對於兩個連續變項的「相關係數」，那剛剛的皮爾森相關其實已經提供了很好的結果
- 並不是只有連續變項能做相關性檢定，二元類別變項或序位變項也可以做相關係數檢定
– 非連續變項的相關係數最常使用的是『Polychoric correlation』
- 很遺憾的，R裡面的內建函數並沒有辦法做這個運算，因此我們需要問問google該怎麼辦。
第五節：相關係數(2)
- 首先，我們必須學會如何安裝套件，使用Rstudio的話，畫面右下角有一個分頁『Packages』，這個分頁可以拿來安裝新的『Package』。在R裡面光發佈在官方網站的就有超過7000個套件，而這些套件都有不一樣的擴充功能，其中套件『ploycor』就是拿來計算Polychoric correlation的套件。
– 我們可以透過下列方法安裝這個套件
– 或是使用函數『install.packages()』安裝指定套件
install.packages("polycor")
- 在安裝完成後，我們必須使用函數『library()』載入這個套件的功能。需要注意的是，當我們執行R的時候，這些額外的套件是不會被載入的，所以每次開啟R以後都必須再次使用函數『library()』(但不需要重新安裝)。
第五節：相關係數(3)
- 安裝完成後，我們需要閱讀套件說明。所有發佈在官方網站的套件都有制式的說明書，我們需要學會如何讀懂他！
第五節：相關係數(4)
– 如果你看不懂每一行的指令在寫甚麼，請一行一行看，並了解該函數的input的格式是什麼。
if(require(mvtnorm)){
set.seed(12345)
data <- rmvnorm(1000, c(0, 0), matrix(c(1, .5, .5, 1), 2, 2))
x <- data[,1]
y <- data[,2]
cor(x, y) # sample correlation
}
## [1] 0.5263698
if(require(mvtnorm)){
x <- cut(x, c(-Inf, .75, Inf))
y <- cut(y, c(-Inf, -1, .5, 1.5, Inf))
polychor(x, y) # 2-step estimate
}
## [1] 0.5230474
if(require(mvtnorm)){
set.seed(12345)
polychor(x, y, ML=TRUE, std.err=TRUE) # ML estimate
}
##
## Polychoric Correlation, ML est. = 0.5231 (0.03819)
## Test of bivariate normality: Chisquare = 2.739, df = 2, p = 0.2543
##
## Row Threshold
## Threshold Std.Err.
## 0.7537 0.04403
##
##
## Column Thresholds
## Threshold Std.Err.
## 1 -0.9842 0.04746
## 2 0.4841 0.04127
## 3 1.5010 0.06118
第五節：相關係數(5)
- 看不懂他的範例沒關係，我們可以多google相關資訊，先了解這個統計方法大概是怎麼算的。
- 簡單來說，Polychoric correlation是假設原來有兩個『連續變項』，兩者服從bivariate normal distribution。
– 但是可能因為臨床意義等原因(如空腹血糖126以上被定義為糖尿病)，使這兩個變項被置換成類別/序位變項。
– 因此Polychoric correlation的想法是，想辦法將兩個類別變項還原成連續變項，再進行相關係數檢定。
– 所以在範例的程式碼，第一個部分是真的製造兩個連續變項，並做Pearson correlation；第二個部分則是強硬的把他切割成數份，並計算Polychoric correlation看答案是不是類似
if(require(mvtnorm)){
set.seed(12345)
data <- rmvnorm(1000, c(0, 0), matrix(c(1, .5, .5, 1), 2, 2))
x <- data[,1]
y <- data[,2]
cor(x, y) # sample correlation
}
## [1] 0.5263698
if(require(mvtnorm)){
x <- cut(x, c(-Inf, .75, Inf))
y <- cut(y, c(-Inf, -1, .5, 1.5, Inf))
polychor(x, y) # 2-step estimate
}
## [1] 0.5230474
第五節：相關係數(6)
- 另外，套件裡的另外一個函數「polyserial()」則是做『類別vs連續』的相關係數，我們就不多做說明，直接套到我們的資料試試：
Result5 = polychor(dat[,"GENDER"], dat[,"LVD"])
Result5
## [1] 0.2317155
Result6 = polyserial(dat[,"K"], dat[,"LVD"])
Result6
## [1] 0.06453313
課程小結
- 本次課程中同學學習到在R語言內進行基本的推論性統計，雖然相關原理並不清楚，但已經學會如何快速的操作R語言進行相對應的分析
– 這節課教的方法裡面大約能面對95%的簡易分析上，之後的統計方法開始會帶有「預測」的功能，我們會將其歸類為「人工智慧」的範疇
– 除此之外還有一些檢定可以用於相依樣本上，像是paired t test及Mcnemar test，你是否能透過google出來呢?
- 這節課最重要的是我們學會了『Package』的使用，這對於我們未來的應用非常重要！