คิดมาก ไม่คิดมาก เรื่องยาก เรื่องง่าย ^[C]

บันทึกนี้ผมนำมาจากเรื่องที่เขียน Note เอาไว้ใน Facebook ครับ เห็นควรว่าน่าจะนำมาแบ่งปันในพื้นที่ชานชลาแชร์ เพราะการเข้าถึง Note ใน FB นั้นค่อนข้างเข้าถึงได้จำกัด เฉพาะกลุ่ม เฉพาะบุคคลเท่านั้น

ในบันทึกนี้มีคำขยายความเพิ่มเติมจากที่เขียนไว้ใน Note นิดหน่อยครับ

เป็นเรื่องที่เกี่ยวกับการพิจารณาข้อมูลครับ

มีข้อมูลอยู่ชุดหนึ่งครับ

ตัวอย่างข้อมูล 10 ค่าแรก

> head(dim1x,10)

 [1] 25 22 23 24 21 23 22 24 24 23

สรุปข้อมูลชุดนี้

> summary(dim1x)   

Min. 1st Qu.  Median    Mean 3rd Qu.    Max.  

19.00   23.00   24.00   23.72   25.00   29.00 

พบว่าข้อมูลชุดนี้อยู่ในช่วง 19.00 (Min.) ถึง 29.00 (Max.) หรือมีช่วงห่างของข้อมูล (พิสัย - Range) เท่ากับ 29 -19 = 10 เท่านั้นเอง ข้อมูลนี้มีค่าอยู่ในช่วงแคบ ๆ เท่านั้น ข้อมูลชุดนี้มีค่าเฉลี่ย (Mean) เท่ากับ 23.72

> length(dim1x)

[1] 439

จำนวนข้อมูลทั้งหมดมีจำนวน 439 ค่า

คราวนี้เรามาดูกันว่า distribution ของข้อมูลชุดนี้จะเป็นแบบไหน หรือย่างไร

ที่เราต้องดู distribution ของข้อมูลก็เพื่อที่จะเลือกใช้วิธีการทางสถิติให้ถูกต้องกับ distribution ของข้อมูล ซึ่งปกติทั่วไป สถิติที่เป็น parametric statistics หรือวิธีการทางสถิติหลายวิธี มีข้อตกลงเบื้องต้นของการทดสอบว่า ข้อมูลที่นำมาทดสอบด้วยวิธีการทางสถิตินี้จะต้องมี distribution ที่เป็นแบบ Normal distribution (Gaussian Distribution)

เช่นค่าเฉลี่ยเป็นค่าการวัดแนวโน้มเข้าสู่ส่วนกลางที่ดีก็ต่อเมื่อ ข้อมูลชุดนั้น ๆ มีการกระจายที่เป็นแบบ Normal distribution

คราวนี้เรามาตรวจสอบ distribution ของข้อมูลชุดข้างต้นว่าเป็นอย่างไร

เริ่มจาก box plots (box and whisker plot) นะครับ

>boxplot(dim1x)

boxplots เป็นวิธีการดู distribution ของข้อมูลที่เป็นแบบ classic วิธีหนึ่งครับ ดูง่ายๆ ว่าข้อมูลเป็น Normal distribution ไหม ก็ดูว่า box plots นั้นสมมาตรหรือไม่ อย่างในรูปข้างบน ถ้าพับครึ่งโดยใช้เส้นทึบ ๆ ในกล่องเป็นจุดกึ่งกลาง ส่วนที่อยู่ข้างบนเส้นทึบกับส่วนที่อยู่ด้านล่างเส้นทึบต้องแนบสนิทกันพอดี

boxplot ข้างต้นจะเห็นข้อมูลที่เป็น outlier อยู่ทั้งหนวด (whisker) เส้นบนและหนวด (whisker) เส้นล่างของ box

คราวนี้มาลอง histogram + probability curve นะครับว่าเป็นอย่างไรบ้าง

>hist(dim1x,prob=T)
>curve(dnorm(x,mean(dim1x),sd(dim1x)),add=T, col="blue",lwd=2)

เขาว่า histogram เฉย ๆ มันขึ้นอยู่กับจำนวนแท่ง (bin) โดยเฉพาะถ้าข้อมูลมี outlier ก็จะทำให้ histogram พาลเบ้ไปซะงั้น แน่ล่ะครับก็ไม่มี one size fits all

แต่ histogram ข้างต้นนั้น เป็น histogram ของ prob (density) เช่นเดียวกับเส้นโค้งระฆังคว่ำนั่นครับ

คิดง่ายๆ density คือ รวมทุกอย่างแล้วได้เท่ากับ 1 หรือคือ probability นั่นเอง

เส้นโค้งปกติหรือที่เรียกว่าเส้นโค้งรูประฆังคว่ำ ถ้าเส้นโค้งนี้เป็นรูประฆังสวย ๆ (ระฆังจะสูงหรือจะเตี้ยก็ได้) ก็ถือว่าการกระจายของข้อมูลเป็นแบบปกติ (แบบโค้งระฆังคว่ำ) แต่ถ้าดูแล้วระฆังมันบุบ ๆ เบี้ยว แหว่ง ๆ เอน ๆ เอียงไปจากปกติ ก็จะถือว่าการกระจายของข้อมูลชุดนี้ไม่เป็นแบบปกติ (จะเป็นแบบไหนนั้นอีกเรื่องนึง)

ถ้าดูจาก histogram + density curve เราจะพบว่า การกระจายของข้อมูลชุดนี้ไม่เบี่ยงเบนไปจากการกระจายแบบ Normal Distribution มากนัก

มาดูที่ plot เขาแนะนำให้ใช้กันครับ

QQ-Norm + QQ-Line

>qqnorm(dim1x, col="blue")
>qqline(dim1x, col="red",lwd=2)

ดูว่าข้อมูลมีการกระจายที่เป็นแบบ normal distribution หรือไม่ก็ให้ดูจุดสีน้ำเงิน (qqnorm - ข้อมูลของเรา - ผมป้ายสีให้เป็นสีน้ำเงิน) ว่าอยู่บน หรือเป็นไปในลักษณะเดียวกับ เส้น (สมมติ - qqline) สีแดง (ผมป้ายให้เป็นสีแดง) หรือไม่

หรือหากข้อมูลมีการกระจายของข้อมูลที่เป็นแบบ normal distribution ก็จะมีการกระจายของข้อมูลเป็นไปในลักษณะของเส้น (สมมติ) สีแดง

คราวนี้เราเห็นชัดขึ้นว่าข้อมูลของเรามีลักษณะเป็น tiles ค่อนข้างมาก

ลองดู ecdf ของข้อมูลครับ (ecdf = Empirical Cumulative Distribution Function)

>plot(ecdf(dim1x),col="blue")

จาก ecdf เราจะเห็นข้อมูลที่มีลักษณะเป็น tile เหมือนกับ qqnorm+ qqline

อ่านเพิ่มเติมเรื่อ ecdf ของข้อมูลได้จาก Wikipedia ครับ หรือจากที่นี่ครับ

หลังจากที่ดูด้วยตา พิจารณาจากใจแล้ว เรามีใช้วิธีการทดสอบทางสถิติดูครับ

test ที่เขานิยมใช้กันมาก และเขาบอกว่ามี power of test มาก ก็คือ Shapiro Wilk test

> shapiro.test(dim1x)

         Shapiro-Wilk normality test

data:  dim1xW = 0.97059, p-value = 1.013e-07

ให้ผลการทดสอบออกมา significant, p-value จากการทดสอบอยู่ที่ประมาณ 1.013/10000000

แต่มีบางคนเขาว่าไว้ว่า  shapiro wilk test นี้จะให้ผลการทดสอบไปในทาง nonsignificant ถ้า n น้อย ๆ (p-value -> 1) แต่ถ้า n เยอะ ๆ หรือมี outlier อยู่ ผลการทดสอบก็จะมีแนวโน้มไปทาง significant (p-value -> 0) หรือ Shapiro Wilk test ค่อนข้างจะเป็น test ที่ หัวโบราณ/ อนุรักษ์นิยม (conservative) สักหน่อย

ลอง test อื่นดูบ้างครับ

> AndersonDarlingTest(dim1x)

      Anderson-Darling test of goodness-of-fit   

     Null hypothesis: uniform distribution

data:  dim1xAn = Inf, p-value = 1.367e-06

test นี้ก็ให้ผลไปในทาง significant

> CramerVonMisesTest(dim1x)

         Cramer-von Mises normality test

data:  dim1xW = 1.172, p-value = 7.37e-10

Warning message:In CramerVonMisesTest(dim1x) :  p-value is smaller than 7.37e-10, cannot be computed more accurately

test นี้ก็ให้ผลไปในทาง significant

> ShapiroFranciaTest(dim1x)

        Shapiro-Francia normality test

data:  dim1xW = 0.97178, p-value = 7.94e-07

test นี้ก็ significant

> JarqueBeraTest(dim1x)

         Robust Jarque Bera Test

 data:  dim1xX-squared = 4.4657, df = 2, p-value = 0.1072

อุ๊บ !!!!! test นี้ให้ผลตรงข้ามกัน กลับเป็น nonsignificant

> ks.test(dim1x, "pnorm", mean(dim1x), sd(dim1x), alternative="two.sided")

         One-sample Kolmogorov-Smirnov test

data:  dim1xD = 0.14718, p-value = 1.098e-08

alternative hypothesis: two-sided

Warning message:In ks.test(dim1x, "pnorm", mean(dim1x), sd(dim1x), alternative = "two.sided") :  ties should not be present for the Kolmogorov-Smirnov test

Smirnov test นี่ เห็นมีใช้ใน SPSS ให้ผลการทดสอบ siginificant เหมือนกัน แต่มีคำเตือนนิดหน่อย

> LillieTest(dim1x)

        Lilliefors (Kolmogorov-Smirnov) normality test

data:  dim1xD = 0.14718, p-value < 2.2e-16

Smirnov test อีก test หนึ่ง test นี้ให้ผลการทดสอบออกมา ค่า p-value -> 0 มากกว่า test อื่น ๆ ที่ผ่านมา

อีก test นึงนะครับ

> PearsonTest(dim1x, n.classes=6, adjust=TRUE)

         Pearson chi-square normality test

data:  dim1xP = 40.221, p-value = 9.566e-09

test นี้ก็ให้ผล significant 

สรุปโดยรวมแล้ว test ต่าง ๆ ให้ผลการทดสอบไปในทาง significant หรือ distribution ของข้อมูลชุดนี้ไม่เเป็นแบบ normal distribution (gaussian distribution) และให้ค่า p-value จากการทดสอบ อยู๋ในช่วงใกล้ ๆ กัน ยกเว้น Lilliefors (Kolmogorov-Smirnov) normality test ที่ให้ค่า p-value -> 0 มากกว่า test อื่นๆ ส่วน  Robust Jarque Bera Test นั้นให้ผลในทิศทางตรงข้าม

เขาอีกนั่นแหละครับ แนะนำว่า 

เอ้ย อย่าไปคิดมากกับมันเลย ไม่ใช่ทุกอย่างที่เราต้องไปดู distribution ของข้อมูล ดูไปดูมาเครียดเปล่า ๆ (ผมว่าเอง ตรงนี้ ฮา)  test หลายอย่างก็มี robust ของ test เองอยู่แล้วด้วย

robust นี่ให้คิดถึงกาแฟเอาไว้นะครับ

robusta, arabica, paracetamol

เข้ม ๆ ขม ๆ ใส่น้ำลงไปหน่อยก็ยังเข้ม ๆ ขม ๆ อยู่ เหมือนบาง test ที่ถึงแม้ข้อมูลมันจะเบี้ยว จะเบ้หรือไม่เป็นไปตาม assumption ที่ตั้งไว้ (เช่นข่อมูลต้องมี distribution ที่เป็นแบบ normal distribution) การทดสอบด้วย test นั้น ๆ ยังคงให้ผลถูกต้องอยู่

อันนี้ต้องไปเรียนเรื่อง robust analyis เอาอีกหัวข้อหนึ่งเพิ่มเติม

ง่าย ๆ ก็ใช้วิธีดูด้วยตา พิจารณาจากใจนั้นแหละครับ ไม่เครียดดี เดี๋ยวจะเข้าทำนอง

"ยักตื้นติดกึก ยักลึกติดกัก"

ฮาาาาาาาาาา

เอ่อนะ ... ท่าจะดี

แล้วสรุปว่าข้อมูลชุดนี้มี distribution แบบไหนอ่ะ ....

อิอิอิ

เราเอง

เพลง: บ่นแก้กลุ้ม
ศิลปิน: สมัย อ่อนวงศ์
คำร้อง: กานท์ การุณวงศ์